source: LMDZ6/trunk/libf/dyn3dmem/advtrac_loc.f90 @ 5813

Last change on this file since 5813 was 5748, checked in by dcugnet, 5 months ago
  • Use REAL(KIND=REAL32) and REAL(KIND=REAL64) Iinstead of REAL and DOUBLE PRECISION

to avoid ambiguity problems in generic procedure when reals are promoted to doubles.

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