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

Last change on this file since 5270 was 5268, checked in by abarral, 2 days ago

.f90 <-> .F90 depending on cpp key use
  • 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.8 KB
RevLine 
[4056]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   !
[4143]10   USE infotrac,     ONLY: nqtot, tracers
[4056]11   USE control_mod,  ONLY: iapp_tracvl, day_step, planet_type
12   USE comconst_mod, ONLY: dtvr
13   USE parallel_lmdz
14   USE Write_Field_loc
15   USE Write_Field
16   USE Bands
17   USE mod_hallo
18   USE Vampir
19   USE times
20   USE advtrac_mod, ONLY: finmasse
[5258]21   USE lmdz_cppkeys_wrapper, ONLY: CPPKEY_DEBUGIO
22   USE strings_mod, ONLY: int2str
[4056]23   IMPLICIT NONE
24   !
25   include "dimensions.h"
26   include "paramet.h"
27   include "comdissip.h"
28   include "comgeom2.h"
29   include "description.h"
30!   include "iniprint.h"
[1632]31
[4056]32   !---------------------------------------------------------------------------
33   !     Arguments
34   !---------------------------------------------------------------------------
35   REAL, INTENT(IN) ::  pbarug(ijb_u:ije_u,llm)
36   REAL, INTENT(IN) ::  pbarvg(ijb_v:ije_v,llm)
37   REAL, INTENT(IN) ::      wg(ijb_u:ije_u,llm)
38   REAL, INTENT(IN) ::       p(ijb_u:ije_u,llmp1)
39   REAL, INTENT(IN) ::  massem(ijb_u:ije_u,llm)
40   REAL, INTENT(INOUT) ::    q(ijb_u:ije_u,llm,nqtot)
41   REAL, INTENT(IN) ::    teta(ijb_u:ije_u,llm)
42   REAL, INTENT(IN) ::      pk(ijb_u:ije_u,llm)
43   !---------------------------------------------------------------------------
44   !     Ajout PPM
45   !---------------------------------------------------------------------------
46   REAL :: massebx(ijb_u:ije_u,llm), masseby(ijb_v:ije_v,llm)
47   !---------------------------------------------------------------------------
48   !     Variables locales
49   !---------------------------------------------------------------------------
[4064]50   INTEGER :: ij, l, iq, iadv
[4056]51   REAL(KIND=KIND(1.d0)) :: t_initial, t_final, tps_cpu
52   REAL :: zdp(ijb_u:ije_u), zdpmin, zdpmax
53   INTEGER, SAVE :: iadvtr=0
54!$OMP THREADPRIVATE(iadvtr)
55   EXTERNAL  minmax
[1632]56
[4056]57   !---------------------------------------------------------------------------
58   !     Rajouts pour PPM
59   !---------------------------------------------------------------------------
60   INTEGER :: indice, n
61   REAL :: dtbon                       ! Pas de temps adaptatif pour que CFL<1
62   REAL :: CFLmaxz, aaa, bbb           ! CFL maximum
63   REAL, DIMENSION(iim,jjb_u:jje_u,llm) :: unatppm, vnatppm, fluxwppm
64   REAL ::    qppm(iim*jjnb_u,llm,nqtot)
65   REAL ::   psppm(iim,jjb_u:jje_u)    ! pression  au sol
66   REAL, DIMENSION(llmp1) :: apppm, bpppm
67   LOGICAL, SAVE :: dum=.TRUE., fill=.TRUE.
68   INTEGER :: ijb, ije, ijbu, ijbv, ijeu, ijev, j
69   TYPE(Request),SAVE :: testRequest
[1848]70!$OMP THREADPRIVATE(testRequest)
[1632]71
[4056]72! Test sur l'eventuelle creation de valeurs negatives de la masse
73   ijb = ij_begin; IF(pole_nord) ijb = ij_begin+iip1
74   ije = ij_end;   IF(pole_sud)  ije = ij_end-iip1
[1632]75
[4056]76!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
77   DO l=1,llm-1
78      DO ij = ijb+1,ije
79         zdp(ij) = pbarug(ij-1,l)    - pbarug(ij,l) &
80                 - pbarvg(ij-iip1,l) + pbarvg(ij,l) &
81                 +     wg(ij,l+1)    -     wg(ij,l)
82      END DO
83! ym  ---> pourquoi jjm-1 et non jjm ? a cause du pole ?
84!     CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 )
85      DO ij = ijb,ije-iip1+1,iip1
86         zdp(ij)=zdp(ij+iip1-1)
87      END DO
88      DO ij = ijb,ije
89         zdp(ij)= zdp(ij)*dtvr/ massem(ij,l)
90      END DO
91!     CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
92! ym ---> eventuellement a revoir
93      CALL minmax( ije-ijb+1, zdp(ijb), zdpmin,zdpmax )
94      IF(MAX(ABS(zdpmin),ABS(zdpmax)) >0.5) &
95         WRITE(*,*)'WARNING DP/P l=',l,'  MIN:',zdpmin,'   MAX:', zdpmax
96   END DO
97!$OMP END DO NOWAIT
[1632]98
[4056]99   !---------------------------------------------------------------------------
100   !   Advection proprement dite (Modification Le Croller (07/2001)
101   !---------------------------------------------------------------------------
[1632]102
[4056]103   !---------------------------------------------------------------------------
104   !   Calcul des moyennes basees sur la masse
105   !---------------------------------------------------------------------------
106!ym   CALL massbar_p(massem,massebx,masseby)
[4058]107!ym   ----> Normalement, inutile pour les schemas classiques
108!ym   ----> Reverifier lors de la parallelisation des autres schemas
[1632]109
[5258]110IF (CPPKEY_DEBUGIO) THEN
[4056]111   CALL WriteField_u('massem',massem)
112   CALL WriteField_u('wg',wg)
113   CALL WriteField_u('pbarug',pbarug)
114   CALL WriteField_v('pbarvg',pbarvg)
115   CALL WriteField_u('p_tmp',p)
116   CALL WriteField_u('pk_tmp',pk)
117   CALL WriteField_u('teta_tmp',teta)
118   DO iq=1,nqtot
119      CALL WriteField_u('q_adv'//trim(int2str(iq)),q(:,:,iq))
120   END DO
[5258]121END IF
[1632]122
123!         
[4056]124!  CALL Register_Hallo_v(pbarvg,llm,1,1,1,1,TestRequest)
125!  CALL SendRequest(TestRequest)
126!!$OMP BARRIER
127!  CALL WaitRequest(TestRequest)
128!$OMP BARRIER
[4058]129
[4056]130!  WRITE(*,*) 'advtrac 157: appel de vlspltgen_loc'
[4058]131   CALL vlspltgen_loc(q, 2., massem, wg, pbarug, pbarvg, dtvr, p, pk, teta )
[1632]132
[5258]133IF (CPPKEY_DEBUGIO) THEN
[4056]134   DO iq = 1, nqtot
135      CALL WriteField_u('q_adv'//trim(int2str(iq)),q(:,:,iq))
136   END DO
[5258]137END IF
[1632]138         
[4056]139   GOTO 1234     
140   !-------------------------------------------------------------------------
141   !       Appel des sous programmes d'advection
142   !-------------------------------------------------------------------------
143   DO iq = 1, nqtot
144!     CALL clock(t_initial)
145      IF(tracers(iq)%parent /= 'air') CYCLE
146      iadv = tracers(iq)%iadv
147      !-----------------------------------------------------------------------
148      SELECT CASE(iadv)
149      !-----------------------------------------------------------------------
150         CASE(0); CYCLE
151         !--------------------------------------------------------------------
152         CASE(10)  !--- Schema de Van Leer I MUSCL
153         !--------------------------------------------------------------------
154!           WRITE(*,*) 'advtrac 239: iq,q(1721,19,:)=',iq,q(1721,19,:)     
155!LF         CALL vlsplt_p(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
[1632]156
[4056]157         !--------------------------------------------------------------------
158         CASE(14)  !--- Schema "pseuDO amont" + test sur humidite specifique
159                   !--- pour la vapeur d'eau. F. Codron
160         !--------------------------------------------------------------------
161!           WRITE(*,*) 'advtrac 248: iq,q(1721,19,:)=',iq,q(1721,19,:)
[4469]162            CALL abort_gcm("advtrac","appel a vlspltqs :schema non parallelise",1)
[4056]163!LF         CALL vlspltqs_p(q(1,1,1),2.,massem,wg,pbarug,pbarvg,dtvr,p,pk,teta )
[1632]164
[4056]165         !--------------------------------------------------------------------
166         CASE(12)  !--- Schema de Frederic Hourdin
167         !--------------------------------------------------------------------
[4469]168            CALL abort_gcm("advtrac","appel a vlspltqs :schema non parallelise",1)
[4056]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
[1632]174
[4056]175         !--------------------------------------------------------------------
176         CASE(13)  !--- Pas de temps adaptatif
177         !--------------------------------------------------------------------
[4469]178            CALL abort_gcm("advtrac","schema non parallelise",1)
[4056]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
[1632]184
[4056]185         !--------------------------------------------------------------------
186         CASE(20)  !--- Schema de pente SLOPES
187         !--------------------------------------------------------------------
[4469]188            CALL abort_gcm("advtrac","schema SLOPES non parallelise",1)
[4056]189            CALL pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
[1632]190
[4056]191         !--------------------------------------------------------------------
192         CASE(30)  !--- Schema de Prather
193         !--------------------------------------------------------------------
[4469]194            CALL abort_gcm("advtrac","schema prather non parallelise",1)
[4056]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)
[1632]199
[4056]200         !--------------------------------------------------------------------
201         CASE(11,16,17,18)   !--- Schemas PPM Lin et Rood
202         !--------------------------------------------------------------------
[4469]203            CALL abort_gcm("advtrac","schema PPM non parallelise",1)
[4056]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)
[1632]224
[4056]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      !----------------------------------------------------------------------
[1632]259
[4056]260      !----------------------------------------------------------------------
[4058]261      ! On impose une seule valeur du traceur au pole Sud j=jjm+1=jjp1 et Nord j=1
[4056]262      !----------------------------------------------------------------------
263      !  CALL traceurpole(q(1,1,iq),massem)
[1632]264
[4056]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
[1632]269
[4056]270   END DO
[1632]271
[4056]2721234 CONTINUE
273!$OMP BARRIER
274   IF(planet_type=="earth") THEN
[1632]275      ijb=ij_begin
276      ije=ij_end
[4056]277!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
278      DO l = 1, llm
[1632]279         DO ij = ijb, ije
[4056]280            finmasse(ij,l) =  p(ij,l) - p(ij,l+1)
281         END DO
282      END DO
283!$OMP END DO
[1632]284
[4056]285      CALL qminimum_loc( q, nqtot, finmasse )
[1632]286
[4056]287   END IF ! of if (planet_type=="earth")
[1632]288
[4056]289END SUBROUTINE advtrac_loc
[1632]290
Note: See TracBrowser for help on using the repository browser.