Ignore:
Timestamp:
Sep 11, 2024, 6:03:07 PM (2 months ago)
Author:
abarral
Message:

Encapsulate files in modules

File:
1 moved

Legend:

Unmodified
Added
Removed

  • LMDZ6/branches/Amaury_dev/libf/dyn3d/lmdz_advtrac.f90

    r5185 r5186  
    1 ! $Id$
    2 
    3 SUBROUTINE advtrac(pbaru, pbarv, p, masse, q, iapptrac, teta, flxw, pk)
    4   !     Auteur :  F. Hourdin
    5 
    6   !     Modif. P. Le Van     (20/12/97)
    7   !            F. Codron     (10/99)
    8   !            D. Le Croller (07/2001)
    9   !            M.A Filiberti (04/2002)
    10 
    11   USE lmdz_infotrac, ONLY: nqtot, tracers, isoCheck
    12   USE control_mod, ONLY: iapp_tracvl, day_step
    13   USE comconst_mod, ONLY: dtvr
    14   USE lmdz_cppkeys_wrapper, ONLY: CPPKEY_DEBUGIO
    15   USE lmdz_strings, ONLY: int2str
    16   USE lmdz_description, ONLY: descript
    17   USE lmdz_libmath, ONLY: minmax
    18   USE lmdz_iniprint, ONLY: lunout, prt_level
    19   USE lmdz_ssum_scopy, ONLY: scopy
    20   USE lmdz_comdissip, ONLY: coefdis, tetavel, tetatemp, gamdissip, niterdis
    21   USE lmdz_comgeom2
    22   USE lmdz_groupe, ONLY: groupe
    23 
    24   USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm
    25   USE lmdz_paramet
    26 IMPLICIT NONE
    27 
    28 
    29 
    30 
    31   !---------------------------------------------------------------------------
    32   !     Arguments
    33   !---------------------------------------------------------------------------
    34   INTEGER, INTENT(OUT) :: iapptrac
    35   REAL, INTENT(IN) :: pbaru(ip1jmp1, llm)
    36   REAL, INTENT(IN) :: pbarv(ip1jm, llm)
    37   REAL, INTENT(INOUT) :: q(ip1jmp1, llm, nqtot)
    38   REAL, INTENT(IN) :: masse(ip1jmp1, llm)
    39   REAL, INTENT(IN) :: p(ip1jmp1, llmp1)
    40   REAL, INTENT(IN) :: teta(ip1jmp1, llm)
    41   REAL, INTENT(IN) :: pk(ip1jmp1, llm)
    42   REAL, INTENT(OUT) :: flxw(ip1jmp1, llm)
    43   !---------------------------------------------------------------------------
    44   !     Ajout PPM
    45   !---------------------------------------------------------------------------
    46   REAL :: massebx(ip1jmp1, llm), masseby(ip1jm, llm)
    47   !---------------------------------------------------------------------------
    48   !     Variables locales
    49   !---------------------------------------------------------------------------
    50   INTEGER :: ij, l, iq, iadv
    51   !   REAL(KIND=KIND(1.d0)) :: t_initial, t_final, tps_cpu
    52   REAL :: zdp(ip1jmp1), zdpmin, zdpmax
    53   INTEGER, SAVE :: iadvtr = 0
    54   REAL, DIMENSION(ip1jmp1, llm) :: pbaruc, pbarug, massem, wg
    55   REAL, DIMENSION(ip1jm, llm) :: pbarvc, pbarvg
    56   SAVE massem, pbaruc, pbarvc
    57   !---------------------------------------------------------------------------
    58   !     Rajouts pour PPM
    59   !---------------------------------------------------------------------------
    60   INTEGER indice, n
    61   REAL :: dtbon                       ! Pas de temps adaptatif pour que CFL<1
    62   REAL :: CFLmaxz, aaa, bbb           ! CFL maximum
    63   REAL, DIMENSION(iim, jjp1, llm) :: unatppm, vnatppm, fluxwppm
    64   REAL :: qppm(iim * jjp1, llm, nqtot)
    65   REAL :: psppm(iim, jjp1)           ! pression  au sol
    66   REAL, DIMENSION(llmp1) :: apppm, bpppm
    67   LOGICAL, SAVE :: dum = .TRUE., fill = .TRUE.
    68 
    69   INTEGER, SAVE :: countcfl = 0
    70   REAL, DIMENSION(ip1jmp1, llm) :: cflx, cflz
    71   REAL, DIMENSION(ip1jm, llm) :: cfly
    72   REAL, DIMENSION(llm), SAVE :: cflxmax, cflymax, cflzmax
    73 
    74   IF(iadvtr == 0) THEN
    75     pbaruc(:, :) = 0
    76     pbarvc(:, :) = 0
    77   END IF
    78 
    79   !--- Accumulation des flux de masse horizontaux
    80   DO l = 1, llm
    81     DO ij = 1, ip1jmp1
    82       pbaruc(ij, l) = pbaruc(ij, l) + pbaru(ij, l)
    83     END DO
    84     DO ij = 1, ip1jm
    85       pbarvc(ij, l) = pbarvc(ij, l) + pbarv(ij, l)
    86     END DO
    87   END DO
    88 
    89   !--- Selection de la masse instantannee des mailles avant le transport.
    90   IF(iadvtr == 0) THEN
    91     CALL SCOPY(ip1jmp1 * llm, masse, 1, massem, 1)
    92     ! CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 )
    93   END IF
    94 
    95   iadvtr = iadvtr + 1
    96   iapptrac = iadvtr
    97 
    98   !--- Test pour savoir si on advecte a ce pas de temps
    99   IF(iadvtr /= iapp_tracvl) RETURN
    100 
    101   !   ..  Modif P.Le Van  ( 20/12/97 )  ....
    102 
    103   !   traitement des flux de masse avant advection.
    104   !       1. calcul de w
    105   !       2. groupement des mailles pres du pole.
    106 
    107   CALL groupe(pbaruc, pbarvc, pbarug, pbarvg, wg)
    108 
    109   !--- Flux de masse diaganostiques traceurs
    110   flxw = wg / REAL(iapp_tracvl)
    111 
    112   !--- Test sur l'eventuelle creation de valeurs negatives de la masse
    113   DO l = 1, llm - 1
    114     DO ij = iip2 + 1, ip1jm
    115       zdp(ij) = pbarug(ij - 1, l) - pbarug(ij, l) &
    116               - pbarvg(ij - iip1, l) + pbarvg(ij, l) &
    117               + wg(ij, l + 1) - wg(ij, l)
    118     END DO
    119     ! ym  ---> pourquoi jjm-1 et non jjm ? a cause du pole ?
    120     CALL SCOPY(jjm - 1, zdp(iip1 + iip1), iip1, zdp(iip2), iip1)
    121     DO ij = iip2, ip1jm
    122       zdp(ij) = zdp(ij) * dtvr / massem(ij, l)
    123     END DO
    124 
    125     CALL minmax (ip1jm - iip1, zdp(iip2), zdpmin, zdpmax)
    126 
    127     IF(MAX(ABS(zdpmin), ABS(zdpmax)) > 0.5) &
    128             WRITE(*, *)'WARNING DP/P l=', l, '  MIN:', zdpmin, ' MAX:', zdpmax
    129 
    130   END DO
    131 
    132   !-------------------------------------------------------------------------
    133   ! Calcul des criteres CFL en X, Y et Z
    134   !-------------------------------------------------------------------------
    135   IF(countcfl == 0.) THEN
    136     cflxmax(:) = 0.
    137     cflymax(:) = 0.
    138     cflzmax(:) = 0.
    139   END IF
    140 
    141   countcfl = countcfl + iapp_tracvl
    142   cflx(:, :) = 0.
    143   cfly(:, :) = 0.
    144   cflz(:, :) = 0.
    145   DO l = 1, llm
    146     DO ij = iip2, ip1jm - 1
    147       IF(pbarug(ij, l)>=0.) THEN
    148         cflx(ij, l) = pbarug(ij, l) * dtvr / masse(ij, l)
    149       ELSE
    150         cflx(ij, l) = -pbarug(ij, l) * dtvr / masse(ij + 1, l)
    151       END IF
    152     END DO
    153   END DO
    154 
    155   DO l = 1, llm
    156     DO ij = iip2, ip1jm - 1, iip1
    157       cflx(ij + iip1, l) = cflx(ij, l)
    158     END DO
    159   END DO
    160 
    161   DO l = 1, llm
    162     DO ij = 1, ip1jm
    163       IF(pbarvg(ij, l)>=0.) THEN
    164         cfly(ij, l) = pbarvg(ij, l) * dtvr / masse(ij, l)
    165       ELSE
    166         cfly(ij, l) = -pbarvg(ij, l) * dtvr / masse(ij + iip1, l)
    167       END IF
    168     END DO
    169   END DO
    170 
    171   DO l = 2, llm
    172     DO ij = 1, ip1jm
    173       IF(wg(ij, l) >= 0.) THEN
    174         cflz(ij, l) = wg(ij, l) * dtvr / masse(ij, l)
    175       ELSE
    176         cflz(ij, l) = -wg(ij, l) * dtvr / masse(ij, l - 1)
    177       END IF
    178     END DO
    179   END DO
    180 
    181   DO l = 1, llm
    182     cflxmax(l) = max(cflxmax(l), maxval(cflx(:, l)))
    183     cflymax(l) = max(cflymax(l), maxval(cfly(:, l)))
    184     cflzmax(l) = max(cflzmax(l), maxval(cflz(:, l)))
    185   END DO
    186 
    187   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    188   ! Par defaut, on sort le diagnostic des CFL tous les jours.
    189   ! Si on veut le sortir a chaque pas d'advection en cas de plantage
    190   !       IF(countcfl==iapp_tracvl) THEN
    191   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    192   IF(countcfl==day_step) THEN
    193     DO l = 1, llm
    194       WRITE(lunout, *) 'L, CFL[xyz]max:', l, cflxmax(l), cflymax(l), cflzmax(l)
    195     END DO
    196     countcfl = 0
    197   END IF
    198 
    199   !---------------------------------------------------------------------------
    200   !   Advection proprement dite (Modification Le Croller (07/2001)
    201   !---------------------------------------------------------------------------
    202 
    203   !---------------------------------------------------------------------------
    204   !   Calcul des moyennes basees sur la masse
    205   !---------------------------------------------------------------------------
    206   CALL massbar(massem, massebx, masseby)
    207 
    208   IF (CPPKEY_DEBUGIO) THEN
    209     CALL WriteField_u('massem', massem)
    210     CALL WriteField_u('wg', wg)
    211     CALL WriteField_u('pbarug', pbarug)
    212     CALL WriteField_v('pbarvg', pbarvg)
    213     CALL WriteField_u('p_tmp', p)
    214     CALL WriteField_u('pk_tmp', pk)
    215     CALL WriteField_u('teta_tmp', teta)
     1MODULE lmdz_advtrac
     2  IMPLICIT NONE; PRIVATE
     3  PUBLIC advtrac
     4
     5CONTAINS
     6
     7  SUBROUTINE advtrac(pbaru, pbarv, p, masse, q, iapptrac, teta, flxw, pk)
     8    !     Auteur :  F. Hourdin
     9
     10    !     Modif. P. Le Van     (20/12/97)
     11    !            F. Codron     (10/99)
     12    !            D. Le Croller (07/2001)
     13    !            M.A Filiberti (04/2002)
     14
     15    USE lmdz_infotrac, ONLY: nqtot, tracers, isoCheck
     16    USE control_mod, ONLY: iapp_tracvl, day_step
     17    USE comconst_mod, ONLY: dtvr
     18    USE lmdz_cppkeys_wrapper, ONLY: CPPKEY_DEBUGIO
     19    USE lmdz_strings, ONLY: int2str
     20    USE lmdz_description, ONLY: descript
     21    USE lmdz_libmath, ONLY: minmax
     22    USE lmdz_iniprint, ONLY: lunout, prt_level
     23    USE lmdz_ssum_scopy, ONLY: scopy
     24    USE lmdz_comdissip, ONLY: coefdis, tetavel, tetatemp, gamdissip, niterdis
     25    USE lmdz_comgeom2
     26    USE lmdz_groupe, ONLY: groupe
     27
     28    USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm
     29    USE lmdz_paramet
     30    USE lmdz_check_isotopes, ONLY: check_isotopes_seq
     31
     32    IMPLICIT NONE
     33
     34
     35
     36
     37    !---------------------------------------------------------------------------
     38    !     Arguments
     39    !---------------------------------------------------------------------------
     40    INTEGER, INTENT(OUT) :: iapptrac
     41    REAL, INTENT(IN) :: pbaru(ip1jmp1, llm)
     42    REAL, INTENT(IN) :: pbarv(ip1jm, llm)
     43    REAL, INTENT(INOUT) :: q(ip1jmp1, llm, nqtot)
     44    REAL, INTENT(IN) :: masse(ip1jmp1, llm)
     45    REAL, INTENT(IN) :: p(ip1jmp1, llmp1)
     46    REAL, INTENT(IN) :: teta(ip1jmp1, llm)
     47    REAL, INTENT(IN) :: pk(ip1jmp1, llm)
     48    REAL, INTENT(OUT) :: flxw(ip1jmp1, llm)
     49    !---------------------------------------------------------------------------
     50    !     Ajout PPM
     51    !---------------------------------------------------------------------------
     52    REAL :: massebx(ip1jmp1, llm), masseby(ip1jm, llm)
     53    !---------------------------------------------------------------------------
     54    !     Variables locales
     55    !---------------------------------------------------------------------------
     56    INTEGER :: ij, l, iq, iadv
     57    !   REAL(KIND=KIND(1.d0)) :: t_initial, t_final, tps_cpu
     58    REAL :: zdp(ip1jmp1), zdpmin, zdpmax
     59    INTEGER, SAVE :: iadvtr = 0
     60    REAL, DIMENSION(ip1jmp1, llm) :: pbaruc, pbarug, massem, wg
     61    REAL, DIMENSION(ip1jm, llm) :: pbarvc, pbarvg
     62    SAVE massem, pbaruc, pbarvc
     63    !---------------------------------------------------------------------------
     64    !     Rajouts pour PPM
     65    !---------------------------------------------------------------------------
     66    INTEGER indice, n
     67    REAL :: dtbon                       ! Pas de temps adaptatif pour que CFL<1
     68    REAL :: CFLmaxz, aaa, bbb           ! CFL maximum
     69    REAL, DIMENSION(iim, jjp1, llm) :: unatppm, vnatppm, fluxwppm
     70    REAL :: qppm(iim * jjp1, llm, nqtot)
     71    REAL :: psppm(iim, jjp1)           ! pression  au sol
     72    REAL, DIMENSION(llmp1) :: apppm, bpppm
     73    LOGICAL, SAVE :: dum = .TRUE., fill = .TRUE.
     74
     75    INTEGER, SAVE :: countcfl = 0
     76    REAL, DIMENSION(ip1jmp1, llm) :: cflx, cflz
     77    REAL, DIMENSION(ip1jm, llm) :: cfly
     78    REAL, DIMENSION(llm), SAVE :: cflxmax, cflymax, cflzmax
     79
     80    IF(iadvtr == 0) THEN
     81      pbaruc(:, :) = 0
     82      pbarvc(:, :) = 0
     83    END IF
     84
     85    !--- Accumulation des flux de masse horizontaux
     86    DO l = 1, llm
     87      DO ij = 1, ip1jmp1
     88        pbaruc(ij, l) = pbaruc(ij, l) + pbaru(ij, l)
     89      END DO
     90      DO ij = 1, ip1jm
     91        pbarvc(ij, l) = pbarvc(ij, l) + pbarv(ij, l)
     92      END DO
     93    END DO
     94
     95    !--- Selection de la masse instantannee des mailles avant le transport.
     96    IF(iadvtr == 0) THEN
     97      CALL SCOPY(ip1jmp1 * llm, masse, 1, massem, 1)
     98      ! CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 )
     99    END IF
     100
     101    iadvtr = iadvtr + 1
     102    iapptrac = iadvtr
     103
     104    !--- Test pour savoir si on advecte a ce pas de temps
     105    IF(iadvtr /= iapp_tracvl) RETURN
     106
     107    !   ..  Modif P.Le Van  ( 20/12/97 )  ....
     108
     109    !   traitement des flux de masse avant advection.
     110    !       1. calcul de w
     111    !       2. groupement des mailles pres du pole.
     112
     113    CALL groupe(pbaruc, pbarvc, pbarug, pbarvg, wg)
     114
     115    !--- Flux de masse diaganostiques traceurs
     116    flxw = wg / REAL(iapp_tracvl)
     117
     118    !--- Test sur l'eventuelle creation de valeurs negatives de la masse
     119    DO l = 1, llm - 1
     120      DO ij = iip2 + 1, ip1jm
     121        zdp(ij) = pbarug(ij - 1, l) - pbarug(ij, l) &
     122                - pbarvg(ij - iip1, l) + pbarvg(ij, l) &
     123                + wg(ij, l + 1) - wg(ij, l)
     124      END DO
     125      ! ym  ---> pourquoi jjm-1 et non jjm ? a cause du pole ?
     126      CALL SCOPY(jjm - 1, zdp(iip1 + iip1), iip1, zdp(iip2), iip1)
     127      DO ij = iip2, ip1jm
     128        zdp(ij) = zdp(ij) * dtvr / massem(ij, l)
     129      END DO
     130
     131      CALL minmax (ip1jm - iip1, zdp(iip2), zdpmin, zdpmax)
     132
     133      IF(MAX(ABS(zdpmin), ABS(zdpmax)) > 0.5) &
     134              WRITE(*, *)'WARNING DP/P l=', l, '  MIN:', zdpmin, ' MAX:', zdpmax
     135
     136    END DO
     137
     138    !-------------------------------------------------------------------------
     139    ! Calcul des criteres CFL en X, Y et Z
     140    !-------------------------------------------------------------------------
     141    IF(countcfl == 0.) THEN
     142      cflxmax(:) = 0.
     143      cflymax(:) = 0.
     144      cflzmax(:) = 0.
     145    END IF
     146
     147    countcfl = countcfl + iapp_tracvl
     148    cflx(:, :) = 0.
     149    cfly(:, :) = 0.
     150    cflz(:, :) = 0.
     151    DO l = 1, llm
     152      DO ij = iip2, ip1jm - 1
     153        IF(pbarug(ij, l)>=0.) THEN
     154          cflx(ij, l) = pbarug(ij, l) * dtvr / masse(ij, l)
     155        ELSE
     156          cflx(ij, l) = -pbarug(ij, l) * dtvr / masse(ij + 1, l)
     157        END IF
     158      END DO
     159    END DO
     160
     161    DO l = 1, llm
     162      DO ij = iip2, ip1jm - 1, iip1
     163        cflx(ij + iip1, l) = cflx(ij, l)
     164      END DO
     165    END DO
     166
     167    DO l = 1, llm
     168      DO ij = 1, ip1jm
     169        IF(pbarvg(ij, l)>=0.) THEN
     170          cfly(ij, l) = pbarvg(ij, l) * dtvr / masse(ij, l)
     171        ELSE
     172          cfly(ij, l) = -pbarvg(ij, l) * dtvr / masse(ij + iip1, l)
     173        END IF
     174      END DO
     175    END DO
     176
     177    DO l = 2, llm
     178      DO ij = 1, ip1jm
     179        IF(wg(ij, l) >= 0.) THEN
     180          cflz(ij, l) = wg(ij, l) * dtvr / masse(ij, l)
     181        ELSE
     182          cflz(ij, l) = -wg(ij, l) * dtvr / masse(ij, l - 1)
     183        END IF
     184      END DO
     185    END DO
     186
     187    DO l = 1, llm
     188      cflxmax(l) = max(cflxmax(l), maxval(cflx(:, l)))
     189      cflymax(l) = max(cflymax(l), maxval(cfly(:, l)))
     190      cflzmax(l) = max(cflzmax(l), maxval(cflz(:, l)))
     191    END DO
     192
     193    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     194    ! Par defaut, on sort le diagnostic des CFL tous les jours.
     195    ! Si on veut le sortir a chaque pas d'advection en cas de plantage
     196    !       IF(countcfl==iapp_tracvl) THEN
     197    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     198    IF(countcfl==day_step) THEN
     199      DO l = 1, llm
     200        WRITE(lunout, *) 'L, CFL[xyz]max:', l, cflxmax(l), cflymax(l), cflzmax(l)
     201      END DO
     202      countcfl = 0
     203    END IF
     204
     205    !---------------------------------------------------------------------------
     206    !   Advection proprement dite (Modification Le Croller (07/2001)
     207    !---------------------------------------------------------------------------
     208
     209    !---------------------------------------------------------------------------
     210    !   Calcul des moyennes basees sur la masse
     211    !---------------------------------------------------------------------------
     212    CALL massbar(massem, massebx, masseby)
     213
     214    IF (CPPKEY_DEBUGIO) THEN
     215      CALL WriteField_u('massem', massem)
     216      CALL WriteField_u('wg', wg)
     217      CALL WriteField_u('pbarug', pbarug)
     218      CALL WriteField_v('pbarvg', pbarvg)
     219      CALL WriteField_u('p_tmp', p)
     220      CALL WriteField_u('pk_tmp', pk)
     221      CALL WriteField_u('teta_tmp', teta)
     222      DO iq = 1, nqtot
     223        CALL WriteField_u('q_adv' // trim(int2str(iq)), q(:, :, iq))
     224      END DO
     225    END IF
     226
     227    IF(isoCheck) WRITE(*, *) 'advtrac 227'
     228    CALL check_isotopes_seq(q, ip1jmp1, 'advtrac 162')
     229
     230    !-------------------------------------------------------------------------
     231    !       Appel des sous programmes d'advection
     232    !-------------------------------------------------------------------------
    216233    DO iq = 1, nqtot
    217       CALL WriteField_u('q_adv' // trim(int2str(iq)), q(:, :, iq))
    218     END DO
    219   END IF
    220 
    221   IF(isoCheck) WRITE(*, *) 'advtrac 227'
    222   CALL check_isotopes_seq(q, ip1jmp1, 'advtrac 162')
    223 
    224   !-------------------------------------------------------------------------
    225   !       Appel des sous programmes d'advection
    226   !-------------------------------------------------------------------------
    227   DO iq = 1, nqtot
    228     !     CALL clock(t_initial)
    229     IF(tracers(iq)%parent /= 'air') CYCLE
    230     iadv = tracers(iq)%iadv
    231     !-----------------------------------------------------------------------
    232     SELECT CASE(iadv)
     234      !     CALL clock(t_initial)
     235      IF(tracers(iq)%parent /= 'air') CYCLE
     236      iadv = tracers(iq)%iadv
    233237      !-----------------------------------------------------------------------
    234     CASE(0); CYCLE
    235     !--------------------------------------------------------------------
    236     CASE(10)  !--- Schema de Van Leer I MUSCL
     238      SELECT CASE(iadv)
     239        !-----------------------------------------------------------------------
     240      CASE(0); CYCLE
    237241      !--------------------------------------------------------------------
    238       !           WRITE(*,*) 'advtrac 239: iq,q(1721,19,:)=',iq,q(1721,19,:)
    239       CALL vlsplt(q, 2., massem, wg, pbarug, pbarvg, dtvr, iq)
    240 
    241       !--------------------------------------------------------------------
    242     CASE(14)  !--- Schema "pseuDO amont" + test sur humidite specifique
    243       !--- pour la vapeur d'eau. F. Codron
    244       !--------------------------------------------------------------------
    245       !           WRITE(*,*) 'advtrac 248: iq,q(1721,19,:)=',iq,q(1721,19,:)
    246       CALL vlspltqs(q, 2., massem, wg, pbarug, pbarvg, dtvr, p, pk, teta, iq)
    247 
    248       !--------------------------------------------------------------------
    249     CASE(12)  !--- Schema de Frederic Hourdin
    250       !--------------------------------------------------------------------
    251       CALL adaptdt(iadv, dtbon, n, pbarug, massem)   ! pas de temps adaptatif
    252       IF(n > 1) WRITE(*, *) 'WARNING horizontal dt=', dtbon, 'dtvr=', dtvr, 'n=', n
    253       DO indice = 1, n
    254         CALL advn(q(1, 1, iq), massem, wg, pbarug, pbarvg, dtbon, 1)
    255       END DO
    256 
    257       !--------------------------------------------------------------------
    258     CASE(13)  !--- Pas de temps adaptatif
    259       !--------------------------------------------------------------------
    260       CALL adaptdt(iadv, dtbon, n, pbarug, massem)
    261       IF(n > 1) WRITE(*, *) 'WARNING horizontal dt=', dtbon, 'dtvr=', dtvr, 'n=', n
    262       DO indice = 1, n
    263         CALL advn(q(1, 1, iq), massem, wg, pbarug, pbarvg, dtbon, 2)
    264       END DO
    265 
    266       !--------------------------------------------------------------------
    267     CASE(20)  !--- Schema de pente SLOPES
    268       !--------------------------------------------------------------------
    269       CALL pentes_ini (q(1, 1, iq), wg, massem, pbarug, pbarvg, 0)
    270 
    271       !--------------------------------------------------------------------
    272     CASE(30)  !--- Schema de Prather
    273       !--------------------------------------------------------------------
    274       ! Pas de temps adaptatif
    275       CALL adaptdt(iadv, dtbon, n, pbarug, massem)
    276       IF(n > 1) WRITE(*, *) 'WARNING horizontal dt=', dtbon, 'dtvr=', dtvr, 'n=', n
    277       CALL prather(q(1, 1, iq), wg, massem, pbarug, pbarvg, n, dtbon)
    278 
    279       !--------------------------------------------------------------------
    280     CASE(11, 16, 17, 18)   !--- Schemas PPM Lin et Rood
    281       !--------------------------------------------------------------------
    282       ! Test sur le flux horizontal
    283       CALL adaptdt(iadv, dtbon, n, pbarug, massem)   ! pas de temps adaptatif
    284       IF(n > 1) WRITE(*, *) 'WARNING horizontal dt=', dtbon, 'dtvr=', dtvr, 'n=', n
    285       ! Test sur le flux vertical
    286       CFLmaxz = 0.
    287       DO l = 2, llm
    288         DO ij = iip2, ip1jm
    289           aaa = wg(ij, l) * dtvr / massem(ij, l)
    290           CFLmaxz = max(CFLmaxz, aaa)
    291           bbb = -wg(ij, l) * dtvr / massem(ij, l - 1)
    292           CFLmaxz = max(CFLmaxz, bbb)
     242      CASE(10)  !--- Schema de Van Leer I MUSCL
     243        !--------------------------------------------------------------------
     244        !           WRITE(*,*) 'advtrac 239: iq,q(1721,19,:)=',iq,q(1721,19,:)
     245        CALL vlsplt(q, 2., massem, wg, pbarug, pbarvg, dtvr, iq)
     246
     247        !--------------------------------------------------------------------
     248      CASE(14)  !--- Schema "pseuDO amont" + test sur humidite specifique
     249        !--- pour la vapeur d'eau. F. Codron
     250        !--------------------------------------------------------------------
     251        !           WRITE(*,*) 'advtrac 248: iq,q(1721,19,:)=',iq,q(1721,19,:)
     252        CALL vlspltqs(q, 2., massem, wg, pbarug, pbarvg, dtvr, p, pk, teta, iq)
     253
     254        !--------------------------------------------------------------------
     255      CASE(12)  !--- Schema de Frederic Hourdin
     256        !--------------------------------------------------------------------
     257        CALL adaptdt(iadv, dtbon, n, pbarug, massem)   ! pas de temps adaptatif
     258        IF(n > 1) WRITE(*, *) 'WARNING horizontal dt=', dtbon, 'dtvr=', dtvr, 'n=', n
     259        DO indice = 1, n
     260          CALL advn(q(1, 1, iq), massem, wg, pbarug, pbarvg, dtbon, 1)
    293261        END DO
    294       END DO
    295       IF(CFLmaxz>=1) WRITE(*, *) 'WARNING vertical', 'CFLmaxz=', CFLmaxz
    296       !----------------------------------------------------------------
    297       !     Ss-prg interface LMDZ.4->PPM3d (ss-prg de Lin)
    298       !----------------------------------------------------------------
    299       CALL interpre(q(1, 1, iq), qppm(1, 1, iq), wg, fluxwppm, massem, &
    300               apppm, bpppm, massebx, masseby, pbarug, pbarvg, &
    301               unatppm, vnatppm, psppm)
    302 
    303       !----------------------------------------------------------------
    304       DO indice = 1, n     !--- VL (version PPM) horiz. et PPM vert.
    305         !----------------------------------------------------------------
    306         SELECT CASE(iadv)
    307           !----------------------------------------------------------
    308         CASE(11)
    309           !----------------------------------------------------------
    310           CALL ppm3d(1, qppm(1, 1, iq), psppm, psppm, unatppm, vnatppm, fluxwppm, dtbon, &
    311                   2, 2, 2, 1, iim, jjp1, 2, llm, apppm, bpppm, 0.01, 6400000, fill, dum, 220.)
    312           !----------------------------------------------------------
    313         CASE(16) !--- Monotonic PPM
    314           !----------------------------------------------------------
    315           CALL ppm3d(1, qppm(1, 1, iq), psppm, psppm, unatppm, vnatppm, fluxwppm, dtbon, &
    316                   3, 3, 3, 1, iim, jjp1, 2, llm, apppm, bpppm, 0.01, 6400000, fill, dum, 220.)
    317           !----------------------------------------------------------
    318         CASE(17) !--- Semi monotonic PPM
    319           !----------------------------------------------------------
    320           CALL ppm3d(1, qppm(1, 1, iq), psppm, psppm, unatppm, vnatppm, fluxwppm, dtbon, &
    321                   4, 4, 4, 1, iim, jjp1, 2, llm, apppm, bpppm, 0.01, 6400000, fill, dum, 220.)
    322           !----------------------------------------------------------
    323         CASE(18) !--- Positive Definite PPM
    324           !----------------------------------------------------------
    325           CALL ppm3d(1, qppm(1, 1, iq), psppm, psppm, unatppm, vnatppm, fluxwppm, dtbon, &
    326                   5, 5, 5, 1, iim, jjp1, 2, llm, apppm, bpppm, 0.01, 6400000, fill, dum, 220.)
    327         END SELECT
    328         !----------------------------------------------------------------
    329       END DO
    330       !----------------------------------------------------------------
    331       !     Ss-prg interface PPM3d-LMDZ.4
    332       !----------------------------------------------------------------
    333       CALL interpost(q(1, 1, iq), qppm(1, 1, iq))
     262
     263        !--------------------------------------------------------------------
     264      CASE(13)  !--- Pas de temps adaptatif
     265        !--------------------------------------------------------------------
     266        CALL adaptdt(iadv, dtbon, n, pbarug, massem)
     267        IF(n > 1) WRITE(*, *) 'WARNING horizontal dt=', dtbon, 'dtvr=', dtvr, 'n=', n
     268        DO indice = 1, n
     269          CALL advn(q(1, 1, iq), massem, wg, pbarug, pbarvg, dtbon, 2)
     270        END DO
     271
     272        !--------------------------------------------------------------------
     273      CASE(20)  !--- Schema de pente SLOPES
     274        !--------------------------------------------------------------------
     275        CALL pentes_ini (q(1, 1, iq), wg, massem, pbarug, pbarvg, 0)
     276
     277        !--------------------------------------------------------------------
     278      CASE(30)  !--- Schema de Prather
     279        !--------------------------------------------------------------------
     280        ! Pas de temps adaptatif
     281        CALL adaptdt(iadv, dtbon, n, pbarug, massem)
     282        IF(n > 1) WRITE(*, *) 'WARNING horizontal dt=', dtbon, 'dtvr=', dtvr, 'n=', n
     283        CALL prather(q(1, 1, iq), wg, massem, pbarug, pbarvg, n, dtbon)
     284
     285        !--------------------------------------------------------------------
     286      CASE(11, 16, 17, 18)   !--- Schemas PPM Lin et Rood
     287        !--------------------------------------------------------------------
     288        ! Test sur le flux horizontal
     289        CALL adaptdt(iadv, dtbon, n, pbarug, massem)   ! pas de temps adaptatif
     290        IF(n > 1) WRITE(*, *) 'WARNING horizontal dt=', dtbon, 'dtvr=', dtvr, 'n=', n
     291        ! Test sur le flux vertical
     292        CFLmaxz = 0.
     293        DO l = 2, llm
     294          DO ij = iip2, ip1jm
     295            aaa = wg(ij, l) * dtvr / massem(ij, l)
     296            CFLmaxz = max(CFLmaxz, aaa)
     297            bbb = -wg(ij, l) * dtvr / massem(ij, l - 1)
     298            CFLmaxz = max(CFLmaxz, bbb)
     299          END DO
     300        END DO
     301        IF(CFLmaxz>=1) WRITE(*, *) 'WARNING vertical', 'CFLmaxz=', CFLmaxz
     302        !----------------------------------------------------------------
     303        !     Ss-prg interface LMDZ.4->PPM3d (ss-prg de Lin)
     304        !----------------------------------------------------------------
     305        CALL interpre(q(1, 1, iq), qppm(1, 1, iq), wg, fluxwppm, massem, &
     306                apppm, bpppm, massebx, masseby, pbarug, pbarvg, &
     307                unatppm, vnatppm, psppm)
     308
     309        !----------------------------------------------------------------
     310        DO indice = 1, n     !--- VL (version PPM) horiz. et PPM vert.
     311          !----------------------------------------------------------------
     312          SELECT CASE(iadv)
     313            !----------------------------------------------------------
     314          CASE(11)
     315            !----------------------------------------------------------
     316            CALL ppm3d(1, qppm(1, 1, iq), psppm, psppm, unatppm, vnatppm, fluxwppm, dtbon, &
     317                    2, 2, 2, 1, iim, jjp1, 2, llm, apppm, bpppm, 0.01, 6400000, fill, dum, 220.)
     318            !----------------------------------------------------------
     319          CASE(16) !--- Monotonic PPM
     320            !----------------------------------------------------------
     321            CALL ppm3d(1, qppm(1, 1, iq), psppm, psppm, unatppm, vnatppm, fluxwppm, dtbon, &
     322                    3, 3, 3, 1, iim, jjp1, 2, llm, apppm, bpppm, 0.01, 6400000, fill, dum, 220.)
     323            !----------------------------------------------------------
     324          CASE(17) !--- Semi monotonic PPM
     325            !----------------------------------------------------------
     326            CALL ppm3d(1, qppm(1, 1, iq), psppm, psppm, unatppm, vnatppm, fluxwppm, dtbon, &
     327                    4, 4, 4, 1, iim, jjp1, 2, llm, apppm, bpppm, 0.01, 6400000, fill, dum, 220.)
     328            !----------------------------------------------------------
     329          CASE(18) !--- Positive Definite PPM
     330            !----------------------------------------------------------
     331            CALL ppm3d(1, qppm(1, 1, iq), psppm, psppm, unatppm, vnatppm, fluxwppm, dtbon, &
     332                    5, 5, 5, 1, iim, jjp1, 2, llm, apppm, bpppm, 0.01, 6400000, fill, dum, 220.)
     333          END SELECT
     334          !----------------------------------------------------------------
     335        END DO
     336        !----------------------------------------------------------------
     337        !     Ss-prg interface PPM3d-LMDZ.4
     338        !----------------------------------------------------------------
     339        CALL interpost(q(1, 1, iq), qppm(1, 1, iq))
     340        !----------------------------------------------------------------------
     341      END SELECT
    334342      !----------------------------------------------------------------------
    335     END SELECT
    336     !----------------------------------------------------------------------
    337 
    338     !----------------------------------------------------------------------
    339     ! On impose une seule valeur du traceur au pole Sud j=jjm+1=jjp1 et Nord j=1
    340     !----------------------------------------------------------------------
    341     !  CALL traceurpole(q(1,1,iq),massem)
    342 
    343     !--- Calcul du temps cpu pour un schema donne
    344     !  CALL clock(t_final)
    345     !ym  tps_cpu=t_final-t_initial
    346     !ym  cpuadv(iq)=cpuadv(iq)+tps_cpu
    347 
    348   END DO
    349 
    350   IF(isoCheck) WRITE(*, *) 'advtrac 402'
    351   CALL check_isotopes_seq(q, ip1jmp1, 'advtrac 397')
    352 
    353   !-------------------------------------------------------------------------
    354   !   on reinitialise a zero les flux de masse cumules
    355   !-------------------------------------------------------------------------
    356   iadvtr = 0
    357 
    358 END SUBROUTINE advtrac
     343
     344      !----------------------------------------------------------------------
     345      ! On impose une seule valeur du traceur au pole Sud j=jjm+1=jjp1 et Nord j=1
     346      !----------------------------------------------------------------------
     347      !  CALL traceurpole(q(1,1,iq),massem)
     348
     349      !--- Calcul du temps cpu pour un schema donne
     350      !  CALL clock(t_final)
     351      !ym  tps_cpu=t_final-t_initial
     352      !ym  cpuadv(iq)=cpuadv(iq)+tps_cpu
     353
     354    END DO
     355
     356    IF(isoCheck) WRITE(*, *) 'advtrac 402'
     357    CALL check_isotopes_seq(q, ip1jmp1, 'advtrac 397')
     358
     359    !-------------------------------------------------------------------------
     360    !   on reinitialise a zero les flux de masse cumules
     361    !-------------------------------------------------------------------------
     362    iadvtr = 0
     363
     364  END SUBROUTINE advtrac
     365
     366END MODULE lmdz_advtrac
Note: See TracChangeset for help on using the changeset viewer.