Ignore:
Timestamp:
Jan 12, 2022, 10:54:09 PM (2 years ago)
Author:
dcugnet
Message:

Most of the changes are intended to help to eventually remove the constraints about the tracers assumptions, in particular water tracers.

  • Remove index tables itr_indice and niadv, replaced by tracers(:)%isAdvected and tracers(:)%isH2OFamily. Most of the loops are now from 1 to nqtot:
    • DO iq=nqo+1,nqtot loops are replaced with: DO iq=1,nqtot

IF(tracers(iq)%isH2Ofamily) CYCLE

  • DO it=1,nbtr; iq=niadv(it+nqo)

and DO it=1,nqtottr; iq=itr_indice(it) loops are replaced with:

it = 0
DO iq = 1, nqtot

IF(.NOT.tracers(iq)%isAdvected .OR. tracers(iq)%isH2Ofamily) CYCLE
it = it+1

  • Move some StratAer? related code from infotrac to infotrac_phy
  • Remove "nqperes" variable:

DO iq=1,nqpere loops are replaced with:
DO iq=1,nqtot

IF(tracers(iq)%parent/='air') CYCLE

  • Cosmetic changes (justification, SELECT CASE instead of multiple IF...) mostly in advtrac* routines.
File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/dyn3d/advtrac.F90

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