Ignore:
Timestamp:
Jul 5, 2011, 10:41:12 AM (13 years ago)
Author:
lguez
Message:

Conversion to free source form, no other change

File:
1 moved

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/dyn3dpar/advtrac_p.F90

    r1492 r1549  
    1 !
    21! $Id$
    3 !
    4 c
    5 c
    6       SUBROUTINE advtrac_p(pbaru,pbarv ,
    7      *                   p,  masse,q,iapptrac,teta,
    8      *                  flxw,
    9      *                  pk   )
    10 
    11 c     Auteur :  F. Hourdin
    12 c
    13 c     Modif. P. Le Van     (20/12/97)
    14 c            F. Codron     (10/99)
    15 c            D. Le Croller (07/2001)
    16 c            M.A Filiberti (04/2002)
    17 c
    18       USE parallel
    19       USE Write_Field_p
    20       USE Bands
    21       USE mod_hallo
    22       USE Vampir
    23       USE times
    24       USE infotrac
    25       USE control_mod
    26       IMPLICIT NONE
    27 c
    28 #include "dimensions.h"
    29 #include "paramet.h"
    30 #include "comconst.h"
    31 #include "comvert.h"
    32 #include "comdissip.h"
    33 #include "comgeom2.h"
    34 #include "logic.h"
    35 #include "temps.h"
    36 #include "ener.h"
    37 #include "description.h"
    38 
    39 c-------------------------------------------------------------------
    40 c     Arguments
    41 c-------------------------------------------------------------------
    42 c     Ajout PPM
    43 c--------------------------------------------------------
    44       REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm)
    45 c--------------------------------------------------------
    46       INTEGER iapptrac
    47       REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
    48       REAL q(ip1jmp1,llm,nqtot),masse(ip1jmp1,llm)
    49       REAL p( ip1jmp1,llmp1 ),teta(ip1jmp1,llm)
    50       REAL pk(ip1jmp1,llm)
    51       REAL               :: flxw(ip1jmp1,llm)
    52 
    53 c-------------------------------------------------------------
    54 c     Variables locales
    55 c-------------------------------------------------------------
    56 
    57       REAL pbaruc(ip1jmp1,llm),pbarvc(ip1jm,llm)
    58       REAL massem(ip1jmp1,llm),zdp(ip1jmp1)
    59       REAL,SAVE::pbarug(ip1jmp1,llm),pbarvg(ip1jm,llm),wg(ip1jmp1,llm)
    60       REAL (kind=kind(1.d0)) :: t_initial, t_final, tps_cpu
    61       INTEGER iadvtr
    62       INTEGER ij,l,iq,iiq
    63       REAL zdpmin, zdpmax
    64       SAVE iadvtr, massem, pbaruc, pbarvc
    65       DATA iadvtr/0/
    66 c$OMP THREADPRIVATE(iadvtr)
    67 c----------------------------------------------------------
    68 c     Rajouts pour PPM
    69 c----------------------------------------------------------
    70       INTEGER indice,n
    71       REAL dtbon ! Pas de temps adaptatif pour que CFL<1
    72       REAL CFLmaxz,aaa,bbb ! CFL maximum
    73       REAL psppm(iim,jjp1) ! pression  au sol
    74       REAL unatppm(iim,jjp1,llm),vnatppm(iim,jjp1,llm)
    75       REAL qppm(iim*jjp1,llm,nqtot)
    76       REAL fluxwppm(iim,jjp1,llm)
    77       REAL apppm(llmp1), bpppm(llmp1)
    78       LOGICAL dum,fill
    79       DATA fill/.true./
    80       DATA dum/.true./
    81       REAL,SAVE :: finmasse(ip1jmp1,llm)
    82       integer ijb,ije,ijb_u,ijb_v,ije_u,ije_v,j
    83       type(Request) :: Request_vanleer
    84       REAL,SAVE :: p_tmp( ip1jmp1,llmp1 )
    85       REAL,SAVE :: teta_tmp(ip1jmp1,llm)
    86       REAL,SAVE :: pk_tmp(ip1jmp1,llm)
    87 
    88       ijb_u=ij_begin
    89       ije_u=ij_end
    90      
    91       ijb_v=ij_begin-iip1
    92       ije_v=ij_end
    93       if (pole_nord) ijb_v=ij_begin
    94       if (pole_sud)  ije_v=ij_end-iip1
    95 
    96       IF(iadvtr.EQ.0) THEN
    97 c         CALL initial0(ijp1llm,pbaruc)
    98 c         CALL initial0(ijmllm,pbarvc)
    99 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)       
    100         DO l=1,llm   
    101           pbaruc(ijb_u:ije_u,l)=0.
    102           pbarvc(ijb_v:ije_v,l)=0.
    103         ENDDO
    104 c$OMP END DO NOWAIT 
    105       ENDIF
    106 
    107 c   accumulation des flux de masse horizontaux
    108 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    109       DO l=1,llm
    110          DO ij = ijb_u,ije_u
    111             pbaruc(ij,l) = pbaruc(ij,l) + pbaru(ij,l)
    112          ENDDO
    113          DO ij = ijb_v,ije_v
    114             pbarvc(ij,l) = pbarvc(ij,l) + pbarv(ij,l)
    115          ENDDO
    116       ENDDO
    117 c$OMP END DO NOWAIT
    118 
    119 c   selection de la masse instantannee des mailles avant le transport.
    120       IF(iadvtr.EQ.0) THEN
    121 
    122 c         CALL SCOPY(ip1jmp1*llm,masse,1,massem,1)
    123           ijb=ij_begin
    124           ije=ij_end
    125 
    126 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    127        DO l=1,llm
    128           massem(ijb:ije,l)=masse(ijb:ije,l)
    129        ENDDO
    130 c$OMP END DO NOWAIT
    131 
    132 ccc         CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 )
    133 c
    134       ENDIF ! of IF(iadvtr.EQ.0)
    135 
    136       iadvtr   = iadvtr+1
    137 
    138 c$OMP MASTER
    139       iapptrac = iadvtr
    140 c$OMP END MASTER
    141 
    142 c   Test pour savoir si on advecte a ce pas de temps
    143 
    144       IF ( iadvtr.EQ.iapp_tracvl ) THEN
    145 c$OMP MASTER
    146         call suspend_timer(timer_caldyn)
    147 c$OMP END MASTER
    148      
    149       ijb=ij_begin
    150       ije=ij_end
    151      
    152 
    153 cc   ..  Modif P.Le Van  ( 20/12/97 )  ....
    154 cc
    155 
    156 c   traitement des flux de masse avant advection.
    157 c     1. calcul de w
    158 c     2. groupement des mailles pres du pole.
    159 
    160         CALL groupe_p( massem, pbaruc,pbarvc, pbarug,pbarvg,wg )
    161 
    162 c$OMP BARRIER
    163 
    164 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    165       DO l=1,llmp1
     2
     3SUBROUTINE advtrac_p(pbaru,pbarv , p,  masse,q,iapptrac,teta, flxw, pk)
     4
     5  !     Auteur :  F. Hourdin
     6  !
     7  !     Modif. P. Le Van     (20/12/97)
     8  !            F. Codron     (10/99)
     9  !            D. Le Croller (07/2001)
     10  !            M.A Filiberti (04/2002)
     11  !
     12  USE parallel
     13  USE Write_Field_p
     14  USE Bands
     15  USE mod_hallo
     16  USE Vampir
     17  USE times
     18  USE infotrac
     19  USE control_mod
     20  IMPLICIT NONE
     21  !
     22  include "dimensions.h"
     23  include "paramet.h"
     24  include "comconst.h"
     25  include "comvert.h"
     26  include "comdissip.h"
     27  include "comgeom2.h"
     28  include "logic.h"
     29  include "temps.h"
     30  include "ener.h"
     31  include "description.h"
     32
     33  !-------------------------------------------------------------------
     34  !     Arguments
     35  !-------------------------------------------------------------------
     36  !     Ajout PPM
     37  !--------------------------------------------------------
     38  REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm)
     39  !--------------------------------------------------------
     40  INTEGER iapptrac
     41  REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
     42  REAL q(ip1jmp1,llm,nqtot),masse(ip1jmp1,llm)
     43  REAL p( ip1jmp1,llmp1 ),teta(ip1jmp1,llm)
     44  REAL pk(ip1jmp1,llm)
     45  REAL               :: flxw(ip1jmp1,llm)
     46
     47  !-------------------------------------------------------------
     48  !     Variables locales
     49  !-------------------------------------------------------------
     50
     51  REAL pbaruc(ip1jmp1,llm),pbarvc(ip1jm,llm)
     52  REAL massem(ip1jmp1,llm),zdp(ip1jmp1)
     53  REAL,SAVE::pbarug(ip1jmp1,llm),pbarvg(ip1jm,llm),wg(ip1jmp1,llm)
     54  REAL (kind=kind(1.d0)) :: t_initial, t_final, tps_cpu
     55  INTEGER iadvtr
     56  INTEGER ij,l,iq,iiq
     57  REAL zdpmin, zdpmax
     58  SAVE iadvtr, massem, pbaruc, pbarvc
     59  DATA iadvtr/0/
     60  !$OMP THREADPRIVATE(iadvtr)
     61  !----------------------------------------------------------
     62  !     Rajouts pour PPM
     63  !----------------------------------------------------------
     64  INTEGER indice,n
     65  REAL dtbon ! Pas de temps adaptatif pour que CFL<1
     66  REAL CFLmaxz,aaa,bbb ! CFL maximum
     67  REAL psppm(iim,jjp1) ! pression  au sol
     68  REAL unatppm(iim,jjp1,llm),vnatppm(iim,jjp1,llm)
     69  REAL qppm(iim*jjp1,llm,nqtot)
     70  REAL fluxwppm(iim,jjp1,llm)
     71  REAL apppm(llmp1), bpppm(llmp1)
     72  LOGICAL dum,fill
     73  DATA fill/.true./
     74  DATA dum/.true./
     75  REAL,SAVE :: finmasse(ip1jmp1,llm)
     76  integer ijb,ije,ijb_u,ijb_v,ije_u,ije_v,j
     77  type(Request) :: Request_vanleer
     78  REAL,SAVE :: p_tmp( ip1jmp1,llmp1 )
     79  REAL,SAVE :: teta_tmp(ip1jmp1,llm)
     80  REAL,SAVE :: pk_tmp(ip1jmp1,llm)
     81
     82  ijb_u=ij_begin
     83  ije_u=ij_end
     84
     85  ijb_v=ij_begin-iip1
     86  ije_v=ij_end
     87  if (pole_nord) ijb_v=ij_begin
     88  if (pole_sud)  ije_v=ij_end-iip1
     89
     90  IF(iadvtr.EQ.0) THEN
     91     !         CALL initial0(ijp1llm,pbaruc)
     92     !         CALL initial0(ijmllm,pbarvc)
     93     !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)       
     94     DO l=1,llm   
     95        pbaruc(ijb_u:ije_u,l)=0.
     96        pbarvc(ijb_v:ije_v,l)=0.
     97     ENDDO
     98     !$OMP END DO NOWAIT 
     99  ENDIF
     100
     101  !   accumulation des flux de masse horizontaux
     102  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
     103  DO l=1,llm
     104     DO ij = ijb_u,ije_u
     105        pbaruc(ij,l) = pbaruc(ij,l) + pbaru(ij,l)
     106     ENDDO
     107     DO ij = ijb_v,ije_v
     108        pbarvc(ij,l) = pbarvc(ij,l) + pbarv(ij,l)
     109     ENDDO
     110  ENDDO
     111  !$OMP END DO NOWAIT
     112
     113  !   selection de la masse instantannee des mailles avant le transport.
     114  IF(iadvtr.EQ.0) THEN
     115
     116     !         CALL SCOPY(ip1jmp1*llm,masse,1,massem,1)
     117     ijb=ij_begin
     118     ije=ij_end
     119
     120     !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
     121     DO l=1,llm
     122        massem(ijb:ije,l)=masse(ijb:ije,l)
     123     ENDDO
     124     !$OMP END DO NOWAIT
     125
     126     !cc         CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 )
     127     !
     128  ENDIF ! of IF(iadvtr.EQ.0)
     129
     130  iadvtr   = iadvtr+1
     131
     132  !$OMP MASTER
     133  iapptrac = iadvtr
     134  !$OMP END MASTER
     135
     136  !   Test pour savoir si on advecte a ce pas de temps
     137
     138  IF ( iadvtr.EQ.iapp_tracvl ) THEN
     139     !$OMP MASTER
     140     call suspend_timer(timer_caldyn)
     141     !$OMP END MASTER
     142
     143     ijb=ij_begin
     144     ije=ij_end
     145
     146
     147     !c   ..  Modif P.Le Van  ( 20/12/97 )  ....
     148     !c
     149
     150     !   traitement des flux de masse avant advection.
     151     !     1. calcul de w
     152     !     2. groupement des mailles pres du pole.
     153
     154     CALL groupe_p( massem, pbaruc,pbarvc, pbarug,pbarvg,wg )
     155
     156     !$OMP BARRIER
     157
     158     !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
     159     DO l=1,llmp1
    166160        p_tmp(ijb:ije,l)=p(ijb:ije,l)
    167       ENDDO
    168 c$OMP END DO NOWAIT
    169      
    170 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    171       DO l=1,llm
     161     ENDDO
     162     !$OMP END DO NOWAIT
     163
     164     !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
     165     DO l=1,llm
    172166        pk_tmp(ijb:ije,l)=pk(ijb:ije,l)
    173167        teta_tmp(ijb:ije,l)=teta(ijb:ije,l)
    174       ENDDO
    175 c$OMP END DO NOWAIT
    176 
    177 c$OMP MASTER
    178       call VTb(VTHallo)
    179 c$OMP END MASTER
    180 
    181       call Register_SwapFieldHallo(pbarug,pbarug,ip1jmp1,llm,
    182      *                             jj_Nb_vanleer,0,0,Request_vanleer)
    183       call Register_SwapFieldHallo(pbarvg,pbarvg,ip1jm,llm,
    184      *                             jj_Nb_vanleer,1,0,Request_vanleer)
    185       call Register_SwapFieldHallo(massem,massem,ip1jmp1,llm,
    186      *                             jj_Nb_vanleer,0,0,Request_vanleer)
    187       call Register_SwapFieldHallo(wg,wg,ip1jmp1,llm,
    188      *                             jj_Nb_vanleer,0,0,Request_vanleer)
    189       call Register_SwapFieldHallo(teta_tmp,teta_tmp,ip1jmp1,llm,
    190      *                             jj_Nb_vanleer,1,1,Request_vanleer)
    191       call Register_SwapFieldHallo(p_tmp,p_tmp,ip1jmp1,llmp1,
    192      *                             jj_Nb_vanleer,1,1,Request_vanleer)
    193       call Register_SwapFieldHallo(pk_tmp,pk_tmp,ip1jmp1,llm,
    194      *                             jj_Nb_vanleer,1,1,Request_vanleer)
    195       do j=1,nqtot
    196         call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm,
    197      *                             jj_nb_vanleer,0,0,Request_vanleer)
    198       enddo
    199 
    200       call SendRequest(Request_vanleer)
    201 c$OMP BARRIER
    202       call WaitRequest(Request_vanleer)
    203 
    204 
    205 c$OMP BARRIER
    206 c$OMP MASTER     
    207       call SetDistrib(jj_nb_vanleer)
    208       call VTe(VTHallo)
    209       call VTb(VTadvection)
    210       call start_timer(timer_vanleer)
    211 c$OMP END MASTER
    212 c$OMP BARRIER
    213      
    214       ! ... Flux de masse diaganostiques traceurs
    215          ijb=ij_begin
    216          ije=ij_end
    217          flxw(ijb:ije,1:llm)=wg(ijb:ije,1:llm)/REAL(iapp_tracvl)
    218 
    219 c  test sur l'eventuelle creation de valeurs negatives de la masse
    220          ijb=ij_begin
    221          ije=ij_end
    222          if (pole_nord) ijb=ij_begin+iip1
    223          if (pole_sud) ije=ij_end-iip1
    224          
    225 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)         
    226          DO l=1,llm-1
    227             DO ij = ijb+1,ije
    228               zdp(ij) =    pbarug(ij-1,l)   - pbarug(ij,l)
    229      s                  - pbarvg(ij-iip1,l) + pbarvg(ij,l)
    230      s                  +       wg(ij,l+1)  - wg(ij,l)
    231             ENDDO
    232            
    233 c            CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 )
    234 c ym  ---> pourquoi jjm-1 et non jjm ? a cause du pole ?
    235            
    236             do ij=ijb,ije-iip1+1,iip1
    237               zdp(ij)=zdp(ij+iip1-1)
    238             enddo
    239            
    240             DO ij = ijb,ije
    241                zdp(ij)= zdp(ij)*dtvr/ massem(ij,l)
    242             ENDDO
    243 
    244 
    245 c            CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
    246 c  ym ---> eventuellement a revoir
    247             CALL minmax ( ije-ijb+1, zdp(ijb), zdpmin,zdpmax )
    248            
    249             IF(MAX(ABS(zdpmin),ABS(zdpmax)).GT.0.5) THEN
    250             PRINT*,'WARNING DP/P l=',l,'  MIN:',zdpmin,
    251      s        '   MAX:', zdpmax
    252             ENDIF
    253 
    254          ENDDO
    255 c$OMP END DO NOWAIT
    256 
    257 c-------------------------------------------------------------------
    258 c   Advection proprement dite (Modification Le Croller (07/2001)
    259 c-------------------------------------------------------------------
    260 
    261 c----------------------------------------------------
    262 c        Calcul des moyennes basées sur la masse
    263 c----------------------------------------------------
    264 
    265 cym      ----> Normalement, inutile pour les schémas classiques
    266 cym      ----> Revérifier lors de la parallélisation des autres schemas
    267    
    268 cym          call massbar_p(massem,massebx,masseby) 
    269 
    270           call vlspltgen_p( q,iadv, 2., massem, wg ,
    271      *                    pbarug,pbarvg,dtvr,p_tmp,pk_tmp,teta_tmp )
    272 
    273          
    274           GOTO 1234     
    275 c-----------------------------------------------------------
    276 c     Appel des sous programmes d'advection
    277 c-----------------------------------------------------------
    278       do iq=1,nqtot
    279 c        call clock(t_initial)
     168     ENDDO
     169     !$OMP END DO NOWAIT
     170
     171     !$OMP MASTER
     172     call VTb(VTHallo)
     173     !$OMP END MASTER
     174
     175     call Register_SwapFieldHallo(pbarug,pbarug,ip1jmp1,llm, &
     176          jj_Nb_vanleer,0,0,Request_vanleer)
     177     call Register_SwapFieldHallo(pbarvg,pbarvg,ip1jm,llm, &
     178          jj_Nb_vanleer,1,0,Request_vanleer)
     179     call Register_SwapFieldHallo(massem,massem,ip1jmp1,llm, &
     180          jj_Nb_vanleer,0,0,Request_vanleer)
     181     call Register_SwapFieldHallo(wg,wg,ip1jmp1,llm, &
     182          jj_Nb_vanleer,0,0,Request_vanleer)
     183     call Register_SwapFieldHallo(teta_tmp,teta_tmp,ip1jmp1,llm, &
     184          jj_Nb_vanleer,1,1,Request_vanleer)
     185     call Register_SwapFieldHallo(p_tmp,p_tmp,ip1jmp1,llmp1, &
     186          jj_Nb_vanleer,1,1,Request_vanleer)
     187     call Register_SwapFieldHallo(pk_tmp,pk_tmp,ip1jmp1,llm, &
     188          jj_Nb_vanleer,1,1,Request_vanleer)
     189     do j=1,nqtot
     190        call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm, &
     191             jj_nb_vanleer,0,0,Request_vanleer)
     192     enddo
     193
     194     call SendRequest(Request_vanleer)
     195     !$OMP BARRIER
     196     call WaitRequest(Request_vanleer)
     197
     198
     199     !$OMP BARRIER
     200     !$OMP MASTER     
     201     call SetDistrib(jj_nb_vanleer)
     202     call VTe(VTHallo)
     203     call VTb(VTadvection)
     204     call start_timer(timer_vanleer)
     205     !$OMP END MASTER
     206     !$OMP BARRIER
     207
     208     ! ... Flux de masse diaganostiques traceurs
     209     ijb=ij_begin
     210     ije=ij_end
     211     flxw(ijb:ije,1:llm)=wg(ijb:ije,1:llm)/REAL(iapp_tracvl)
     212
     213     !  test sur l'eventuelle creation de valeurs negatives de la masse
     214     ijb=ij_begin
     215     ije=ij_end
     216     if (pole_nord) ijb=ij_begin+iip1
     217     if (pole_sud) ije=ij_end-iip1
     218
     219     !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)         
     220     DO l=1,llm-1
     221        DO ij = ijb+1,ije
     222           zdp(ij) =    pbarug(ij-1,l)   - pbarug(ij,l) &
     223                - pbarvg(ij-iip1,l) + pbarvg(ij,l) &
     224                +       wg(ij,l+1)  - wg(ij,l)
     225        ENDDO
     226
     227        !            CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 )
     228        ! ym  ---> pourquoi jjm-1 et non jjm ? a cause du pole ?
     229
     230        do ij=ijb,ije-iip1+1,iip1
     231           zdp(ij)=zdp(ij+iip1-1)
     232        enddo
     233
     234        DO ij = ijb,ije
     235           zdp(ij)= zdp(ij)*dtvr/ massem(ij,l)
     236        ENDDO
     237
     238
     239        !            CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
     240        !  ym ---> eventuellement a revoir
     241        CALL minmax ( ije-ijb+1, zdp(ijb), zdpmin,zdpmax )
     242
     243        IF(MAX(ABS(zdpmin),ABS(zdpmax)).GT.0.5) THEN
     244           PRINT*,'WARNING DP/P l=',l,'  MIN:',zdpmin, &
     245                '   MAX:', zdpmax
     246        ENDIF
     247
     248     ENDDO
     249     !$OMP END DO NOWAIT
     250
     251     !-------------------------------------------------------------------
     252     !   Advection proprement dite (Modification Le Croller (07/2001)
     253     !-------------------------------------------------------------------
     254
     255     !----------------------------------------------------
     256     !        Calcul des moyennes basées sur la masse
     257     !----------------------------------------------------
     258
     259     !ym      ----> Normalement, inutile pour les schémas classiques
     260     !ym      ----> Revérifier lors de la parallélisation des autres schemas
     261
     262     !ym          call massbar_p(massem,massebx,masseby) 
     263
     264     call vlspltgen_p( q,iadv, 2., massem, wg , &
     265          pbarug,pbarvg,dtvr,p_tmp,pk_tmp,teta_tmp )
     266
     267
     268     GOTO 1234     
     269     !-----------------------------------------------------------
     270     !     Appel des sous programmes d'advection
     271     !-----------------------------------------------------------
     272     do iq=1,nqtot
     273        !        call clock(t_initial)
    280274        if(iadv(iq) == 0) cycle
    281 c   ----------------------------------------------------------------
    282 c   Schema de Van Leer I MUSCL
    283 c   ----------------------------------------------------------------
     275        !   ----------------------------------------------------------------
     276        !   Schema de Van Leer I MUSCL
     277        !   ----------------------------------------------------------------
    284278        if(iadv(iq).eq.10) THEN
    285      
    286             call vlsplt_p(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
    287 
    288 c   ----------------------------------------------------------------
    289 c   Schema "pseudo amont" + test sur humidite specifique
    290 C    pour la vapeur d'eau. F. Codron
    291 c   ----------------------------------------------------------------
     279
     280           call vlsplt_p(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
     281
     282           !   ----------------------------------------------------------------
     283           !   Schema "pseudo amont" + test sur humidite specifique
     284           !    pour la vapeur d'eau. F. Codron
     285           !   ----------------------------------------------------------------
    292286        else if(iadv(iq).eq.14) then
    293 c
    294 cym           stop 'advtrac : appel à vlspltqs :schema non parallelise'
    295            CALL vlspltqs_p( q(1,1,1), 2., massem, wg ,
    296      *                 pbarug,pbarvg,dtvr,p_tmp,pk_tmp,teta_tmp )
    297 c   ----------------------------------------------------------------
    298 c   Schema de Frederic Hourdin
    299 c   ----------------------------------------------------------------
     287           !
     288           !ym        stop 'advtrac : appel à vlspltqs :schema non parallelise'
     289           CALL vlspltqs_p( q(1,1,1), 2., massem, wg , &
     290                pbarug,pbarvg,dtvr,p_tmp,pk_tmp,teta_tmp )
     291           !   ----------------------------------------------------------------
     292           !   Schema de Frederic Hourdin
     293           !   ----------------------------------------------------------------
    300294        else if(iadv(iq).eq.12) then
    301           stop 'advtrac : schema non parallelise'
    302 c            Pas de temps adaptatif
     295           stop 'advtrac : schema non parallelise'
     296           !            Pas de temps adaptatif
    303297           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
    304298           if (n.GT.1) then
    305            write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
    306      s             dtvr,'n=',n
     299              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
     300                   dtvr,'n=',n
    307301           endif
    308302           do indice=1,n
    309             call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
     303              call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
    310304           end do
    311305        else if(iadv(iq).eq.13) then
    312           stop 'advtrac : schema non parallelise'
    313 c            Pas de temps adaptatif
     306           stop 'advtrac : schema non parallelise'
     307           !            Pas de temps adaptatif
    314308           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
    315309           if (n.GT.1) then
    316            write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
    317      s             dtvr,'n=',n
     310              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
     311                   dtvr,'n=',n
    318312           endif
    319           do indice=1,n
    320             call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
    321           end do
    322 c   ----------------------------------------------------------------
    323 c   Schema de pente SLOPES
    324 c   ----------------------------------------------------------------
     313           do indice=1,n
     314              call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
     315           end do
     316           !   ----------------------------------------------------------------
     317           !   Schema de pente SLOPES
     318           !   ----------------------------------------------------------------
    325319        else if (iadv(iq).eq.20) then
    326           stop 'advtrac : schema non parallelise'
    327 
    328             call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
    329 
    330 c   ----------------------------------------------------------------
    331 c   Schema de Prather
    332 c   ----------------------------------------------------------------
     320           stop 'advtrac : schema non parallelise'
     321
     322           call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
     323
     324           !   ----------------------------------------------------------------
     325           !   Schema de Prather
     326           !   ----------------------------------------------------------------
    333327        else if (iadv(iq).eq.30) then
    334           stop 'advtrac : schema non parallelise'
    335 c            Pas de temps adaptatif
     328           stop 'advtrac : schema non parallelise'
     329           !            Pas de temps adaptatif
    336330           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
    337331           if (n.GT.1) then
    338            write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
    339      s             dtvr,'n=',n
     332              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
     333                   dtvr,'n=',n
    340334           endif
    341            call  prather(q(1,1,iq),wg,massem,pbarug,pbarvg,
    342      s                     n,dtbon)
    343 c   ----------------------------------------------------------------
    344 c   Schemas PPM Lin et Rood
    345 c   ----------------------------------------------------------------
    346        else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND.
    347      s                     iadv(iq).LE.18)) then
     335           call  prather(q(1,1,iq),wg,massem,pbarug,pbarvg, &
     336                n,dtbon)
     337           !   ----------------------------------------------------------------
     338           !   Schemas PPM Lin et Rood
     339           !   ----------------------------------------------------------------
     340        else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND. &
     341             iadv(iq).LE.18)) then
    348342
    349343           stop 'advtrac : schema non parallelise'
    350344
    351 c        Test sur le flux horizontal
    352 c        Pas de temps adaptatif
    353          call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
    354          if (n.GT.1) then
    355            write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
    356      s             dtvr,'n=',n
    357          endif
    358 c        Test sur le flux vertical
    359          CFLmaxz=0.
    360          do l=2,llm
    361            do ij=iip2,ip1jm
    362             aaa=wg(ij,l)*dtvr/massem(ij,l)
    363             CFLmaxz=max(CFLmaxz,aaa)
    364             bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
    365             CFLmaxz=max(CFLmaxz,bbb)
     345           !        Test sur le flux horizontal
     346           !        Pas de temps adaptatif
     347           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
     348           if (n.GT.1) then
     349              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
     350                   dtvr,'n=',n
     351           endif
     352           !        Test sur le flux vertical
     353           CFLmaxz=0.
     354           do l=2,llm
     355              do ij=iip2,ip1jm
     356                 aaa=wg(ij,l)*dtvr/massem(ij,l)
     357                 CFLmaxz=max(CFLmaxz,aaa)
     358                 bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
     359                 CFLmaxz=max(CFLmaxz,bbb)
     360              enddo
    366361           enddo
    367          enddo
    368          if (CFLmaxz.GE.1) then
    369             write(*,*) 'WARNING vertical','CFLmaxz=', CFLmaxz
    370          endif
    371 
    372 c-----------------------------------------------------------
    373 c        Ss-prg interface LMDZ.4->PPM3d
    374 c-----------------------------------------------------------
    375 
    376           call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem,
    377      s                 apppm,bpppm,massebx,masseby,pbarug,pbarvg,
    378      s                 unatppm,vnatppm,psppm)
    379 
    380           do indice=1,n
    381 c---------------------------------------------------------------------
    382 c                         VL (version PPM) horiz. et PPM vert.
    383 c---------------------------------------------------------------------
    384                 if (iadv(iq).eq.11) then
    385 c                  Ss-prg PPM3d de Lin
    386                   call ppm3d(1,qppm(1,1,iq),
    387      s                       psppm,psppm,
    388      s                       unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1,
    389      s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
    390      s                       fill,dum,220.)
    391 
    392 c----------------------------------------------------------------------
    393 c                           Monotonic PPM
    394 c----------------------------------------------------------------------
    395                else if (iadv(iq).eq.16) then
    396 c                  Ss-prg PPM3d de Lin
    397                   call ppm3d(1,qppm(1,1,iq),
    398      s                       psppm,psppm,
    399      s                       unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1,
    400      s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
    401      s                       fill,dum,220.)
    402 c---------------------------------------------------------------------
    403 
    404 c---------------------------------------------------------------------
    405 c                           Semi Monotonic PPM
    406 c---------------------------------------------------------------------
    407                else if (iadv(iq).eq.17) then
    408 c                  Ss-prg PPM3d de Lin
    409                   call ppm3d(1,qppm(1,1,iq),
    410      s                       psppm,psppm,
    411      s                       unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1,
    412      s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
    413      s                       fill,dum,220.)
    414 c---------------------------------------------------------------------
    415 
    416 c---------------------------------------------------------------------
    417 c                         Positive Definite PPM
    418 c---------------------------------------------------------------------
    419                 else if (iadv(iq).eq.18) then
    420 c                  Ss-prg PPM3d de Lin
    421                   call ppm3d(1,qppm(1,1,iq),
    422      s                       psppm,psppm,
    423      s                       unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1,
    424      s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
    425      s                       fill,dum,220.)
    426 c---------------------------------------------------------------------
    427                 endif
    428             enddo
    429 c-----------------------------------------------------------------
    430 c               Ss-prg interface PPM3d-LMDZ.4
    431 c-----------------------------------------------------------------
    432                   call interpost(q(1,1,iq),qppm(1,1,iq))
    433             endif
    434 c----------------------------------------------------------------------
    435 
    436 c-----------------------------------------------------------------
    437 c On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
    438 c et Nord j=1
    439 c-----------------------------------------------------------------
    440 
    441 c                  call traceurpole(q(1,1,iq),massem)
    442 
    443 c calcul du temps cpu pour un schema donne
    444 
    445 c                  call clock(t_final)
    446 cym                  tps_cpu=t_final-t_initial
    447 cym                  cpuadv(iq)=cpuadv(iq)+tps_cpu
    448 
    449        end DO
    450 
    451 1234  CONTINUE
    452 c$OMP BARRIER
    453 
    454       if (planet_type=="earth") then
     362           if (CFLmaxz.GE.1) then
     363              write(*,*) 'WARNING vertical','CFLmaxz=', CFLmaxz
     364           endif
     365
     366           !-----------------------------------------------------------
     367           !        Ss-prg interface LMDZ.4->PPM3d
     368           !-----------------------------------------------------------
     369
     370           call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem, &
     371                apppm,bpppm,massebx,masseby,pbarug,pbarvg, &
     372                unatppm,vnatppm,psppm)
     373
     374           do indice=1,n
     375              !----------------------------------------------------------------
     376              !                         VL (version PPM) horiz. et PPM vert.
     377              !----------------------------------------------------------------
     378              if (iadv(iq).eq.11) then
     379                 !                  Ss-prg PPM3d de Lin
     380                 call ppm3d(1,qppm(1,1,iq), &
     381                      psppm,psppm, &
     382                      unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1, &
     383                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
     384                      fill,dum,220.)
     385
     386                 !-------------------------------------------------------------
     387                 !                           Monotonic PPM
     388                 !-------------------------------------------------------------
     389              else if (iadv(iq).eq.16) then
     390                 !                  Ss-prg PPM3d de Lin
     391                 call ppm3d(1,qppm(1,1,iq), &
     392                      psppm,psppm, &
     393                      unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1, &
     394                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
     395                      fill,dum,220.)
     396                 !-------------------------------------------------------------
     397
     398                 !-------------------------------------------------------------
     399                 !                           Semi Monotonic PPM
     400                 !-------------------------------------------------------------
     401              else if (iadv(iq).eq.17) then
     402                 !                  Ss-prg PPM3d de Lin
     403                 call ppm3d(1,qppm(1,1,iq), &
     404                      psppm,psppm, &
     405                      unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1, &
     406                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
     407                      fill,dum,220.)
     408                 !-------------------------------------------------------------
     409
     410                 !-------------------------------------------------------------
     411                 !                         Positive Definite PPM
     412                 !-------------------------------------------------------------
     413              else if (iadv(iq).eq.18) then
     414                 !                  Ss-prg PPM3d de Lin
     415                 call ppm3d(1,qppm(1,1,iq), &
     416                      psppm,psppm, &
     417                      unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1, &
     418                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
     419                      fill,dum,220.)
     420                 !-------------------------------------------------------------
     421              endif
     422           enddo
     423           !-----------------------------------------------------------------
     424           !               Ss-prg interface PPM3d-LMDZ.4
     425           !-----------------------------------------------------------------
     426           call interpost(q(1,1,iq),qppm(1,1,iq))
     427        endif
     428        !----------------------------------------------------------------------
     429
     430        !-----------------------------------------------------------------
     431        ! On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
     432        ! et Nord j=1
     433        !-----------------------------------------------------------------
     434
     435        !                  call traceurpole(q(1,1,iq),massem)
     436
     437        ! calcul du temps cpu pour un schema donne
     438
     439        !                  call clock(t_final)
     440        !ym                  tps_cpu=t_final-t_initial
     441        !ym                  cpuadv(iq)=cpuadv(iq)+tps_cpu
     442
     443     end DO
     444
     4451234 CONTINUE
     446     !$OMP BARRIER
     447
     448     if (planet_type=="earth") then
    455449
    456450        ijb=ij_begin
    457451        ije=ij_end
    458452
    459 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
     453        !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
    460454        DO l = 1, llm
    461          DO ij = ijb, ije
    462            finmasse(ij,l) =  p(ij,l) - p(ij,l+1)
    463          ENDDO
     455           DO ij = ijb, ije
     456              finmasse(ij,l) =  p(ij,l) - p(ij,l+1)
     457           ENDDO
    464458        ENDDO
    465 c$OMP END DO
     459        !$OMP END DO
    466460
    467461        CALL qminimum_p( q, 2, finmasse )
    468462
    469 c------------------------------------------------------------------
    470 c   on reinitialise a zero les flux de masse cumules
    471 c---------------------------------------------------
    472 c          iadvtr=0
    473 
    474 c$OMP MASTER
     463        !------------------------------------------------------------------
     464        !   on reinitialise a zero les flux de masse cumules
     465        !---------------------------------------------------
     466        !          iadvtr=0
     467
     468        !$OMP MASTER
    475469        call VTe(VTadvection)
    476470        call stop_timer(timer_vanleer)
    477471        call VTb(VThallo)
    478 c$OMP END MASTER
     472        !$OMP END MASTER
    479473
    480474        do j=1,nqtot
    481           call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm,
    482      *                             jj_nb_caldyn,0,0,Request_vanleer)
     475           call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm, &
     476                jj_nb_caldyn,0,0,Request_vanleer)
    483477        enddo
    484478
    485         call Register_SwapFieldHallo(flxw,flxw,ip1jmp1,llm,
    486      *       jj_nb_caldyn,0,0,Request_vanleer)
     479        call Register_SwapFieldHallo(flxw,flxw,ip1jmp1,llm, &
     480             jj_nb_caldyn,0,0,Request_vanleer)
    487481
    488482        call SendRequest(Request_vanleer)
    489 c$OMP BARRIER
     483        !$OMP BARRIER
    490484        call WaitRequest(Request_vanleer)     
    491485
    492 c$OMP BARRIER
    493 c$OMP MASTER
     486        !$OMP BARRIER
     487        !$OMP MASTER
    494488        call SetDistrib(jj_nb_caldyn)
    495489        call VTe(VThallo)
    496490        call resume_timer(timer_caldyn)
    497 c$OMP END MASTER
    498 c$OMP BARRIER   
    499           iadvtr=0
    500       endif ! of if (planet_type=="earth")
    501        ENDIF ! if iadvtr.EQ.iapp_tracvl
    502 
    503        RETURN
    504        END
    505 
     491 !$OMP END MASTER
     492 !$OMP BARRIER 
     493        iadvtr=0
     494     endif ! of if (planet_type=="earth")
     495  ENDIF ! if iadvtr.EQ.iapp_tracvl
     496
     497END SUBROUTINE advtrac_p
Note: See TracChangeset for help on using the changeset viewer.