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

Last change on this file since 5410 was 5324, checked in by abarral, 5 weeks ago

[WIP] Remove uses of DEBUGIO cpp key (deprecated)
  • 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
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   !
[5281]10   USE comgeom2_mod_h
[5280]11   USE comdissip_mod_h
[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
[5258]23   USE strings_mod, ONLY: int2str
[5271]24   USE dimensions_mod, ONLY: iim, jjm, llm, ndm
[5285]25USE paramet_mod_h
[5271]26IMPLICIT NONE
[4056]27   !
[5271]28
[5272]29
[5282]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
110!         
[4056]111!  CALL Register_Hallo_v(pbarvg,llm,1,1,1,1,TestRequest)
112!  CALL SendRequest(TestRequest)
113!!$OMP BARRIER
114!  CALL WaitRequest(TestRequest)
115!$OMP BARRIER
[4058]116
[4056]117!  WRITE(*,*) 'advtrac 157: appel de vlspltgen_loc'
[4058]118   CALL vlspltgen_loc(q, 2., massem, wg, pbarug, pbarvg, dtvr, p, pk, teta )
[1632]119         
[4056]120   GOTO 1234     
121   !-------------------------------------------------------------------------
122   !       Appel des sous programmes d'advection
123   !-------------------------------------------------------------------------
124   DO iq = 1, nqtot
125!     CALL clock(t_initial)
126      IF(tracers(iq)%parent /= 'air') CYCLE
127      iadv = tracers(iq)%iadv
128      !-----------------------------------------------------------------------
129      SELECT CASE(iadv)
130      !-----------------------------------------------------------------------
131         CASE(0); CYCLE
132         !--------------------------------------------------------------------
133         CASE(10)  !--- Schema de Van Leer I MUSCL
134         !--------------------------------------------------------------------
135!           WRITE(*,*) 'advtrac 239: iq,q(1721,19,:)=',iq,q(1721,19,:)     
136!LF         CALL vlsplt_p(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
[1632]137
[4056]138         !--------------------------------------------------------------------
139         CASE(14)  !--- Schema "pseuDO amont" + test sur humidite specifique
140                   !--- pour la vapeur d'eau. F. Codron
141         !--------------------------------------------------------------------
142!           WRITE(*,*) 'advtrac 248: iq,q(1721,19,:)=',iq,q(1721,19,:)
[4469]143            CALL abort_gcm("advtrac","appel a vlspltqs :schema non parallelise",1)
[4056]144!LF         CALL vlspltqs_p(q(1,1,1),2.,massem,wg,pbarug,pbarvg,dtvr,p,pk,teta )
[1632]145
[4056]146         !--------------------------------------------------------------------
147         CASE(12)  !--- Schema de Frederic Hourdin
148         !--------------------------------------------------------------------
[4469]149            CALL abort_gcm("advtrac","appel a vlspltqs :schema non parallelise",1)
[4056]150            CALL adaptdt(iadv,dtbon,n,pbarug,massem)   ! pas de temps adaptatif
151            IF(n > 1) WRITE(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',dtvr,'n=',n
152            DO indice=1,n
153              CALL advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
154            END DO
[1632]155
[4056]156         !--------------------------------------------------------------------
157         CASE(13)  !--- Pas de temps adaptatif
158         !--------------------------------------------------------------------
[4469]159            CALL abort_gcm("advtrac","schema non parallelise",1)
[4056]160            CALL adaptdt(iadv,dtbon,n,pbarug,massem)
161            IF(n > 1) WRITE(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',dtvr,'n=',n
162            DO indice=1,n
163               CALL advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
164            END DO
[1632]165
[4056]166         !--------------------------------------------------------------------
167         CASE(20)  !--- Schema de pente SLOPES
168         !--------------------------------------------------------------------
[4469]169            CALL abort_gcm("advtrac","schema SLOPES non parallelise",1)
[4056]170            CALL pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
[1632]171
[4056]172         !--------------------------------------------------------------------
173         CASE(30)  !--- Schema de Prather
174         !--------------------------------------------------------------------
[4469]175            CALL abort_gcm("advtrac","schema prather non parallelise",1)
[4056]176            ! Pas de temps adaptatif
177            CALL adaptdt(iadv,dtbon,n,pbarug,massem)
178            IF(n > 1) WRITE(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',dtvr,'n=',n
179            CALL prather(q(1,1,iq),wg,massem,pbarug,pbarvg,n,dtbon)
[1632]180
[4056]181         !--------------------------------------------------------------------
182         CASE(11,16,17,18)   !--- Schemas PPM Lin et Rood
183         !--------------------------------------------------------------------
[4469]184            CALL abort_gcm("advtrac","schema PPM non parallelise",1)
[4056]185            ! Test sur le flux horizontal
186            CALL adaptdt(iadv,dtbon,n,pbarug,massem)   ! pas de temps adaptatif
187            IF(n > 1) WRITE(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',dtvr,'n=',n
188            ! Test sur le flux vertical
189            CFLmaxz=0.
190            DO l=2,llm
191               DO ij=iip2,ip1jm
192                  aaa=wg(ij,l)*dtvr/massem(ij,l)
193                  CFLmaxz=max(CFLmaxz,aaa)
194                  bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
195                  CFLmaxz=max(CFLmaxz,bbb)
196               END DO
197            END DO
198            IF(CFLmaxz.GE.1) WRITE(*,*) 'WARNING vertical','CFLmaxz=', CFLmaxz
199            !----------------------------------------------------------------
200            !     Ss-prg interface LMDZ.4->PPM3d (ss-prg de Lin)
201            !----------------------------------------------------------------
202            CALL interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem, &
203                 apppm,bpppm,massebx,masseby,pbarug,pbarvg, &
204                 unatppm,vnatppm,psppm)
[1632]205
[4056]206            !----------------------------------------------------------------
207            DO indice=1,n     !--- VL (version PPM) horiz. et PPM vert.
208            !----------------------------------------------------------------
209               SELECT CASE(iadv)
210                  !----------------------------------------------------------
211                  CASE(11)
212                  !----------------------------------------------------------
213                     CALL ppm3d(1,qppm(1,1,iq),psppm,psppm,unatppm,vnatppm,fluxwppm,dtbon, &
214                                2,2,2,1,iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,fill,dum,220.)
215                  !----------------------------------------------------------
216                  CASE(16) !--- Monotonic PPM
217                  !----------------------------------------------------------
218                     CALL ppm3d(1,qppm(1,1,iq),psppm,psppm,unatppm,vnatppm,fluxwppm,dtbon, &
219                                3,3,3,1,iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,fill,dum,220.)
220                  !----------------------------------------------------------
221                  CASE(17) !--- Semi monotonic PPM
222                  !----------------------------------------------------------
223                     CALL ppm3d(1,qppm(1,1,iq),psppm,psppm,unatppm,vnatppm,fluxwppm,dtbon, &
224                                4,4,4,1,iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, fill,dum,220.)
225                  !----------------------------------------------------------
226                  CASE(18) !--- Positive Definite PPM
227                  !----------------------------------------------------------
228                     CALL ppm3d(1,qppm(1,1,iq),psppm,psppm,unatppm,vnatppm,fluxwppm,dtbon, &
229                                5,5,5,1,iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,fill,dum,220.)
230               END SELECT
231            !----------------------------------------------------------------
232            END DO
233            !----------------------------------------------------------------
234            !     Ss-prg interface PPM3d-LMDZ.4
235            !----------------------------------------------------------------
236            CALL interpost(q(1,1,iq),qppm(1,1,iq))
237      !----------------------------------------------------------------------
238      END SELECT
239      !----------------------------------------------------------------------
[1632]240
[4056]241      !----------------------------------------------------------------------
[4058]242      ! On impose une seule valeur du traceur au pole Sud j=jjm+1=jjp1 et Nord j=1
[4056]243      !----------------------------------------------------------------------
244      !  CALL traceurpole(q(1,1,iq),massem)
[1632]245
[4056]246      !--- Calcul du temps cpu pour un schema donne
247      !  CALL clock(t_final)
248      !ym  tps_cpu=t_final-t_initial
249      !ym  cpuadv(iq)=cpuadv(iq)+tps_cpu
[1632]250
[4056]251   END DO
[1632]252
[4056]2531234 CONTINUE
254!$OMP BARRIER
255   IF(planet_type=="earth") THEN
[1632]256      ijb=ij_begin
257      ije=ij_end
[4056]258!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
259      DO l = 1, llm
[1632]260         DO ij = ijb, ije
[4056]261            finmasse(ij,l) =  p(ij,l) - p(ij,l+1)
262         END DO
263      END DO
264!$OMP END DO
[1632]265
[4056]266      CALL qminimum_loc( q, nqtot, finmasse )
[1632]267
[4056]268   END IF ! of if (planet_type=="earth")
[1632]269
[4056]270END SUBROUTINE advtrac_loc
[1632]271
Note: See TracBrowser for help on using the repository browser.