source: LMDZ6/trunk/libf/dyn3dmem/advtrac_loc.F90 @ 5182

Last change on this file since 5182 was 4469, checked in by Laurent Fairhead, 20 months ago

Replaced STOP instructions by calls to abort_gcm for a cleaner exit
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
File size: 12.7 KB
RevLine 
[4056]1
[1632]2#define DEBUG_IO
3#undef DEBUG_IO
[4056]4SUBROUTINE advtrac_loc(pbarug, pbarvg, wg, p, massem, q, teta, pk)
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   !
[4143]12   USE infotrac,     ONLY: nqtot, tracers
[4056]13   USE control_mod,  ONLY: iapp_tracvl, day_step, planet_type
14   USE comconst_mod, ONLY: dtvr
15   USE parallel_lmdz
16   USE Write_Field_loc
17   USE Write_Field
18   USE Bands
19   USE mod_hallo
20   USE Vampir
21   USE times
22   USE advtrac_mod, ONLY: finmasse
[1632]23
[4056]24   IMPLICIT NONE
25   !
26   include "dimensions.h"
27   include "paramet.h"
28   include "comdissip.h"
29   include "comgeom2.h"
30   include "description.h"
31!   include "iniprint.h"
[1632]32
[4056]33   !---------------------------------------------------------------------------
34   !     Arguments
35   !---------------------------------------------------------------------------
36   REAL, INTENT(IN) ::  pbarug(ijb_u:ije_u,llm)
37   REAL, INTENT(IN) ::  pbarvg(ijb_v:ije_v,llm)
38   REAL, INTENT(IN) ::      wg(ijb_u:ije_u,llm)
39   REAL, INTENT(IN) ::       p(ijb_u:ije_u,llmp1)
40   REAL, INTENT(IN) ::  massem(ijb_u:ije_u,llm)
41   REAL, INTENT(INOUT) ::    q(ijb_u:ije_u,llm,nqtot)
42   REAL, INTENT(IN) ::    teta(ijb_u:ije_u,llm)
43   REAL, INTENT(IN) ::      pk(ijb_u:ije_u,llm)
44   !---------------------------------------------------------------------------
45   !     Ajout PPM
46   !---------------------------------------------------------------------------
47   REAL :: massebx(ijb_u:ije_u,llm), masseby(ijb_v:ije_v,llm)
48   !---------------------------------------------------------------------------
49   !     Variables locales
50   !---------------------------------------------------------------------------
[4064]51   INTEGER :: ij, l, iq, iadv
[4056]52   REAL(KIND=KIND(1.d0)) :: t_initial, t_final, tps_cpu
53   REAL :: zdp(ijb_u:ije_u), zdpmin, zdpmax
54   INTEGER, SAVE :: iadvtr=0
55!$OMP THREADPRIVATE(iadvtr)
56   EXTERNAL  minmax
[1632]57
[4056]58   !---------------------------------------------------------------------------
59   !     Rajouts pour PPM
60   !---------------------------------------------------------------------------
61   INTEGER :: indice, n
62   REAL :: dtbon                       ! Pas de temps adaptatif pour que CFL<1
63   REAL :: CFLmaxz, aaa, bbb           ! CFL maximum
64   REAL, DIMENSION(iim,jjb_u:jje_u,llm) :: unatppm, vnatppm, fluxwppm
65   REAL ::    qppm(iim*jjnb_u,llm,nqtot)
66   REAL ::   psppm(iim,jjb_u:jje_u)    ! pression  au sol
67   REAL, DIMENSION(llmp1) :: apppm, bpppm
68   LOGICAL, SAVE :: dum=.TRUE., fill=.TRUE.
69   INTEGER :: ijb, ije, ijbu, ijbv, ijeu, ijev, j
70   TYPE(Request),SAVE :: testRequest
[1848]71!$OMP THREADPRIVATE(testRequest)
[1632]72
[4056]73! Test sur l'eventuelle creation de valeurs negatives de la masse
74   ijb = ij_begin; IF(pole_nord) ijb = ij_begin+iip1
75   ije = ij_end;   IF(pole_sud)  ije = ij_end-iip1
[1632]76
[4056]77!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
78   DO l=1,llm-1
79      DO ij = ijb+1,ije
80         zdp(ij) = pbarug(ij-1,l)    - pbarug(ij,l) &
81                 - pbarvg(ij-iip1,l) + pbarvg(ij,l) &
82                 +     wg(ij,l+1)    -     wg(ij,l)
83      END DO
84! ym  ---> pourquoi jjm-1 et non jjm ? a cause du pole ?
85!     CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 )
86      DO ij = ijb,ije-iip1+1,iip1
87         zdp(ij)=zdp(ij+iip1-1)
88      END DO
89      DO ij = ijb,ije
90         zdp(ij)= zdp(ij)*dtvr/ massem(ij,l)
91      END DO
92!     CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
93! ym ---> eventuellement a revoir
94      CALL minmax( ije-ijb+1, zdp(ijb), zdpmin,zdpmax )
95      IF(MAX(ABS(zdpmin),ABS(zdpmax)) >0.5) &
96         WRITE(*,*)'WARNING DP/P l=',l,'  MIN:',zdpmin,'   MAX:', zdpmax
97   END DO
98!$OMP END DO NOWAIT
[1632]99
[4056]100   !---------------------------------------------------------------------------
101   !   Advection proprement dite (Modification Le Croller (07/2001)
102   !---------------------------------------------------------------------------
[1632]103
[4056]104   !---------------------------------------------------------------------------
105   !   Calcul des moyennes basees sur la masse
106   !---------------------------------------------------------------------------
107!ym   CALL massbar_p(massem,massebx,masseby)
[4058]108!ym   ----> Normalement, inutile pour les schemas classiques
109!ym   ----> Reverifier lors de la parallelisation des autres schemas
[1632]110
111#ifdef DEBUG_IO   
[4056]112   CALL WriteField_u('massem',massem)
113   CALL WriteField_u('wg',wg)
114   CALL WriteField_u('pbarug',pbarug)
115   CALL WriteField_v('pbarvg',pbarvg)
116   CALL WriteField_u('p_tmp',p)
117   CALL WriteField_u('pk_tmp',pk)
118   CALL WriteField_u('teta_tmp',teta)
119   DO iq=1,nqtot
120      CALL WriteField_u('q_adv'//trim(int2str(iq)),q(:,:,iq))
121   END DO
[1632]122#endif
123
124!         
[4056]125!  CALL Register_Hallo_v(pbarvg,llm,1,1,1,1,TestRequest)
126!  CALL SendRequest(TestRequest)
127!!$OMP BARRIER
128!  CALL WaitRequest(TestRequest)
129!$OMP BARRIER
[4058]130
[4056]131!  WRITE(*,*) 'advtrac 157: appel de vlspltgen_loc'
[4058]132   CALL vlspltgen_loc(q, 2., massem, wg, pbarug, pbarvg, dtvr, p, pk, teta )
[1632]133
134#ifdef DEBUG_IO     
[4056]135   DO iq = 1, nqtot
136      CALL WriteField_u('q_adv'//trim(int2str(iq)),q(:,:,iq))
137   END DO
[1632]138#endif
139         
[4056]140   GOTO 1234     
141   !-------------------------------------------------------------------------
142   !       Appel des sous programmes d'advection
143   !-------------------------------------------------------------------------
144   DO iq = 1, nqtot
145!     CALL clock(t_initial)
146      IF(tracers(iq)%parent /= 'air') CYCLE
147      iadv = tracers(iq)%iadv
148      !-----------------------------------------------------------------------
149      SELECT CASE(iadv)
150      !-----------------------------------------------------------------------
151         CASE(0); CYCLE
152         !--------------------------------------------------------------------
153         CASE(10)  !--- Schema de Van Leer I MUSCL
154         !--------------------------------------------------------------------
155!           WRITE(*,*) 'advtrac 239: iq,q(1721,19,:)=',iq,q(1721,19,:)     
156!LF         CALL vlsplt_p(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
[1632]157
[4056]158         !--------------------------------------------------------------------
159         CASE(14)  !--- Schema "pseuDO amont" + test sur humidite specifique
160                   !--- pour la vapeur d'eau. F. Codron
161         !--------------------------------------------------------------------
162!           WRITE(*,*) 'advtrac 248: iq,q(1721,19,:)=',iq,q(1721,19,:)
[4469]163            CALL abort_gcm("advtrac","appel a vlspltqs :schema non parallelise",1)
[4056]164!LF         CALL vlspltqs_p(q(1,1,1),2.,massem,wg,pbarug,pbarvg,dtvr,p,pk,teta )
[1632]165
[4056]166         !--------------------------------------------------------------------
167         CASE(12)  !--- Schema de Frederic Hourdin
168         !--------------------------------------------------------------------
[4469]169            CALL abort_gcm("advtrac","appel a vlspltqs :schema non parallelise",1)
[4056]170            CALL adaptdt(iadv,dtbon,n,pbarug,massem)   ! pas de temps adaptatif
171            IF(n > 1) WRITE(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',dtvr,'n=',n
172            DO indice=1,n
173              CALL advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
174            END DO
[1632]175
[4056]176         !--------------------------------------------------------------------
177         CASE(13)  !--- Pas de temps adaptatif
178         !--------------------------------------------------------------------
[4469]179            CALL abort_gcm("advtrac","schema non parallelise",1)
[4056]180            CALL adaptdt(iadv,dtbon,n,pbarug,massem)
181            IF(n > 1) WRITE(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',dtvr,'n=',n
182            DO indice=1,n
183               CALL advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
184            END DO
[1632]185
[4056]186         !--------------------------------------------------------------------
187         CASE(20)  !--- Schema de pente SLOPES
188         !--------------------------------------------------------------------
[4469]189            CALL abort_gcm("advtrac","schema SLOPES non parallelise",1)
[4056]190            CALL pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
[1632]191
[4056]192         !--------------------------------------------------------------------
193         CASE(30)  !--- Schema de Prather
194         !--------------------------------------------------------------------
[4469]195            CALL abort_gcm("advtrac","schema prather non parallelise",1)
[4056]196            ! Pas de temps adaptatif
197            CALL adaptdt(iadv,dtbon,n,pbarug,massem)
198            IF(n > 1) WRITE(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',dtvr,'n=',n
199            CALL prather(q(1,1,iq),wg,massem,pbarug,pbarvg,n,dtbon)
[1632]200
[4056]201         !--------------------------------------------------------------------
202         CASE(11,16,17,18)   !--- Schemas PPM Lin et Rood
203         !--------------------------------------------------------------------
[4469]204            CALL abort_gcm("advtrac","schema PPM non parallelise",1)
[4056]205            ! Test sur le flux horizontal
206            CALL adaptdt(iadv,dtbon,n,pbarug,massem)   ! pas de temps adaptatif
207            IF(n > 1) WRITE(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',dtvr,'n=',n
208            ! Test sur le flux vertical
209            CFLmaxz=0.
210            DO l=2,llm
211               DO ij=iip2,ip1jm
212                  aaa=wg(ij,l)*dtvr/massem(ij,l)
213                  CFLmaxz=max(CFLmaxz,aaa)
214                  bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
215                  CFLmaxz=max(CFLmaxz,bbb)
216               END DO
217            END DO
218            IF(CFLmaxz.GE.1) WRITE(*,*) 'WARNING vertical','CFLmaxz=', CFLmaxz
219            !----------------------------------------------------------------
220            !     Ss-prg interface LMDZ.4->PPM3d (ss-prg de Lin)
221            !----------------------------------------------------------------
222            CALL interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem, &
223                 apppm,bpppm,massebx,masseby,pbarug,pbarvg, &
224                 unatppm,vnatppm,psppm)
[1632]225
[4056]226            !----------------------------------------------------------------
227            DO indice=1,n     !--- VL (version PPM) horiz. et PPM vert.
228            !----------------------------------------------------------------
229               SELECT CASE(iadv)
230                  !----------------------------------------------------------
231                  CASE(11)
232                  !----------------------------------------------------------
233                     CALL ppm3d(1,qppm(1,1,iq),psppm,psppm,unatppm,vnatppm,fluxwppm,dtbon, &
234                                2,2,2,1,iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,fill,dum,220.)
235                  !----------------------------------------------------------
236                  CASE(16) !--- Monotonic PPM
237                  !----------------------------------------------------------
238                     CALL ppm3d(1,qppm(1,1,iq),psppm,psppm,unatppm,vnatppm,fluxwppm,dtbon, &
239                                3,3,3,1,iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,fill,dum,220.)
240                  !----------------------------------------------------------
241                  CASE(17) !--- Semi monotonic PPM
242                  !----------------------------------------------------------
243                     CALL ppm3d(1,qppm(1,1,iq),psppm,psppm,unatppm,vnatppm,fluxwppm,dtbon, &
244                                4,4,4,1,iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, fill,dum,220.)
245                  !----------------------------------------------------------
246                  CASE(18) !--- Positive Definite PPM
247                  !----------------------------------------------------------
248                     CALL ppm3d(1,qppm(1,1,iq),psppm,psppm,unatppm,vnatppm,fluxwppm,dtbon, &
249                                5,5,5,1,iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,fill,dum,220.)
250               END SELECT
251            !----------------------------------------------------------------
252            END DO
253            !----------------------------------------------------------------
254            !     Ss-prg interface PPM3d-LMDZ.4
255            !----------------------------------------------------------------
256            CALL interpost(q(1,1,iq),qppm(1,1,iq))
257      !----------------------------------------------------------------------
258      END SELECT
259      !----------------------------------------------------------------------
[1632]260
[4056]261      !----------------------------------------------------------------------
[4058]262      ! On impose une seule valeur du traceur au pole Sud j=jjm+1=jjp1 et Nord j=1
[4056]263      !----------------------------------------------------------------------
264      !  CALL traceurpole(q(1,1,iq),massem)
[1632]265
[4056]266      !--- Calcul du temps cpu pour un schema donne
267      !  CALL clock(t_final)
268      !ym  tps_cpu=t_final-t_initial
269      !ym  cpuadv(iq)=cpuadv(iq)+tps_cpu
[1632]270
[4056]271   END DO
[1632]272
[4056]2731234 CONTINUE
274!$OMP BARRIER
275   IF(planet_type=="earth") THEN
[1632]276      ijb=ij_begin
277      ije=ij_end
[4056]278!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
279      DO l = 1, llm
[1632]280         DO ij = ijb, ije
[4056]281            finmasse(ij,l) =  p(ij,l) - p(ij,l+1)
282         END DO
283      END DO
284!$OMP END DO
[1632]285
[4056]286      CALL qminimum_loc( q, nqtot, finmasse )
[1632]287
[4056]288   END IF ! of if (planet_type=="earth")
[1632]289
[4056]290END SUBROUTINE advtrac_loc
[1632]291
Note: See TracBrowser for help on using the repository browser.