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

Last change on this file since 4057 was 4056, checked in by dcugnet, 2 years ago

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.
  • 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.6 KB
Line 
1
2#define DEBUG_IO
3#undef DEBUG_IO
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   !
12   USE infotrac,     ONLY: nqtot, tracers,ok_iso_verif
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
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"
32
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   !---------------------------------------------------------------------------
51   INTEGER :: ij, l, iq, iiq, iadv
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
57
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
71!$OMP THREADPRIVATE(testRequest)
72
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
76
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
99
100   !---------------------------------------------------------------------------
101   !   Advection proprement dite (Modification Le Croller (07/2001)
102   !---------------------------------------------------------------------------
103
104   !---------------------------------------------------------------------------
105   !   Calcul des moyennes basees sur la masse
106   !---------------------------------------------------------------------------
107!ym   CALL massbar_p(massem,massebx,masseby)
108!ym   ----> Normalement, inutile pour les schémas classiques
109!ym   ----> Revérifier lors de la parallélisation des autres schemas
110
111#ifdef DEBUG_IO   
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
122#endif
123
124!         
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
130                 
131!  WRITE(*,*) 'advtrac 157: appel de vlspltgen_loc'
132   CALL vlspltgen_loc(q, tracers(:)%iadv, 2., massem, wg, pbarug, pbarvg, dtvr, p, pk, teta )
133
134#ifdef DEBUG_IO     
135   DO iq = 1, nqtot
136      CALL WriteField_u('q_adv'//trim(int2str(iq)),q(:,:,iq))
137   END DO
138#endif
139         
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)
157
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,:)
163            STOP 'advtrac : appel a vlspltqs :schema non parallelise'
164!LF         CALL vlspltqs_p(q(1,1,1),2.,massem,wg,pbarug,pbarvg,dtvr,p,pk,teta )
165
166         !--------------------------------------------------------------------
167         CASE(12)  !--- Schema de Frederic Hourdin
168         !--------------------------------------------------------------------
169            CALL adaptdt(iadv,dtbon,n,pbarug,massem)   ! pas de temps adaptatif
170            IF(n > 1) WRITE(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',dtvr,'n=',n
171            DO indice=1,n
172              CALL advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
173            END DO
174
175         !--------------------------------------------------------------------
176         CASE(13)  !--- Pas de temps adaptatif
177         !--------------------------------------------------------------------
178            STOP 'advtrac : schema non parallelise'
179            CALL adaptdt(iadv,dtbon,n,pbarug,massem)
180            IF(n > 1) WRITE(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',dtvr,'n=',n
181            DO indice=1,n
182               CALL advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
183            END DO
184
185         !--------------------------------------------------------------------
186         CASE(20)  !--- Schema de pente SLOPES
187         !--------------------------------------------------------------------
188            STOP 'advtrac : schema non parallelise'
189            CALL pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
190
191         !--------------------------------------------------------------------
192         CASE(30)  !--- Schema de Prather
193         !--------------------------------------------------------------------
194            STOP 'advtrac : schema non parallelise'
195            ! Pas de temps adaptatif
196            CALL adaptdt(iadv,dtbon,n,pbarug,massem)
197            IF(n > 1) WRITE(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',dtvr,'n=',n
198            CALL prather(q(1,1,iq),wg,massem,pbarug,pbarvg,n,dtbon)
199
200         !--------------------------------------------------------------------
201         CASE(11,16,17,18)   !--- Schemas PPM Lin et Rood
202         !--------------------------------------------------------------------
203            STOP 'advtrac : schema non parallelise'
204            ! Test sur le flux horizontal
205            CALL adaptdt(iadv,dtbon,n,pbarug,massem)   ! pas de temps adaptatif
206            IF(n > 1) WRITE(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',dtvr,'n=',n
207            ! Test sur le flux vertical
208            CFLmaxz=0.
209            DO l=2,llm
210               DO ij=iip2,ip1jm
211                  aaa=wg(ij,l)*dtvr/massem(ij,l)
212                  CFLmaxz=max(CFLmaxz,aaa)
213                  bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
214                  CFLmaxz=max(CFLmaxz,bbb)
215               END DO
216            END DO
217            IF(CFLmaxz.GE.1) WRITE(*,*) 'WARNING vertical','CFLmaxz=', CFLmaxz
218            !----------------------------------------------------------------
219            !     Ss-prg interface LMDZ.4->PPM3d (ss-prg de Lin)
220            !----------------------------------------------------------------
221            CALL interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem, &
222                 apppm,bpppm,massebx,masseby,pbarug,pbarvg, &
223                 unatppm,vnatppm,psppm)
224
225            !----------------------------------------------------------------
226            DO indice=1,n     !--- VL (version PPM) horiz. et PPM vert.
227            !----------------------------------------------------------------
228               SELECT CASE(iadv)
229                  !----------------------------------------------------------
230                  CASE(11)
231                  !----------------------------------------------------------
232                     CALL ppm3d(1,qppm(1,1,iq),psppm,psppm,unatppm,vnatppm,fluxwppm,dtbon, &
233                                2,2,2,1,iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,fill,dum,220.)
234                  !----------------------------------------------------------
235                  CASE(16) !--- Monotonic PPM
236                  !----------------------------------------------------------
237                     CALL ppm3d(1,qppm(1,1,iq),psppm,psppm,unatppm,vnatppm,fluxwppm,dtbon, &
238                                3,3,3,1,iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,fill,dum,220.)
239                  !----------------------------------------------------------
240                  CASE(17) !--- Semi monotonic PPM
241                  !----------------------------------------------------------
242                     CALL ppm3d(1,qppm(1,1,iq),psppm,psppm,unatppm,vnatppm,fluxwppm,dtbon, &
243                                4,4,4,1,iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, fill,dum,220.)
244                  !----------------------------------------------------------
245                  CASE(18) !--- Positive Definite PPM
246                  !----------------------------------------------------------
247                     CALL ppm3d(1,qppm(1,1,iq),psppm,psppm,unatppm,vnatppm,fluxwppm,dtbon, &
248                                5,5,5,1,iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,fill,dum,220.)
249               END SELECT
250            !----------------------------------------------------------------
251            END DO
252            !----------------------------------------------------------------
253            !     Ss-prg interface PPM3d-LMDZ.4
254            !----------------------------------------------------------------
255            CALL interpost(q(1,1,iq),qppm(1,1,iq))
256      !----------------------------------------------------------------------
257      END SELECT
258      !----------------------------------------------------------------------
259
260      !----------------------------------------------------------------------
261      ! On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1 et Nord j=1
262      !----------------------------------------------------------------------
263      !  CALL traceurpole(q(1,1,iq),massem)
264
265      !--- Calcul du temps cpu pour un schema donne
266      !  CALL clock(t_final)
267      !ym  tps_cpu=t_final-t_initial
268      !ym  cpuadv(iq)=cpuadv(iq)+tps_cpu
269
270   END DO
271
2721234 CONTINUE
273!$OMP BARRIER
274   IF(planet_type=="earth") THEN
275      ijb=ij_begin
276      ije=ij_end
277!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
278      DO l = 1, llm
279         DO ij = ijb, ije
280            finmasse(ij,l) =  p(ij,l) - p(ij,l+1)
281         END DO
282      END DO
283!$OMP END DO
284
285      CALL qminimum_loc( q, nqtot, finmasse )
286
287   END IF ! of if (planet_type=="earth")
288
289END SUBROUTINE advtrac_loc
290
Note: See TracBrowser for help on using the repository browser.