source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/config/ppsrc/dyn/vlsplt_p.f @ 224

Last change on this file since 224 was 224, checked in by ymipsl, 10 years ago
File size: 56.0 KB
Line 
1      SUBROUTINE vlsplt_p(q,pente_max,masse,w,pbaru,pbarv,pdt)
2c
3c     Auteurs:   P.Le Van, F.Hourdin, F.Forget 
4c
5c    ********************************************************************
6c     Shema  d'advection " pseudo amont " .
7c    ********************************************************************
8c     q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....
9c
10c   pente_max facteur de limitation des pentes: 2 en general
11c                                               0 pour un schema amont
12c   pbaru,pbarv,w flux de masse en u ,v ,w
13c   pdt pas de temps
14c
15c   --------------------------------------------------------------------
16      USE parallel_lmdz
17      USE mod_hallo
18      USE Vampir
19      IMPLICIT NONE
20c
21!-----------------------------------------------------------------------
22!   INCLUDE 'dimensions.h'
23!
24!   dimensions.h contient les dimensions du modele
25!   ndm est tel que iim=2**ndm
26!-----------------------------------------------------------------------
27
28      INTEGER iim,jjm,llm,ndm
29
30      PARAMETER (iim= 128,jjm=96,llm=64,ndm=1)
31
32!-----------------------------------------------------------------------
33!
34! $Header$
35!
36!
37!  ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
38!                 veillez  n'utiliser que des ! pour les commentaires
39!                 et  bien positionner les & des lignes de continuation
40!                 (les placer en colonne 6 et en colonne 73)
41!
42!
43!-----------------------------------------------------------------------
44!   INCLUDE 'paramet.h'
45
46      INTEGER  iip1,iip2,iip3,jjp1,llmp1,llmp2,llmm1
47      INTEGER  kftd,ip1jm,ip1jmp1,ip1jmi1,ijp1llm
48      INTEGER  ijmllm,mvar
49      INTEGER jcfil,jcfllm
50
51      PARAMETER( iip1= iim+1,iip2=iim+2,iip3=iim+3                       &
52     &    ,jjp1=jjm+1-1/jjm)
53      PARAMETER( llmp1 = llm+1,  llmp2 = llm+2, llmm1 = llm-1 )
54      PARAMETER( kftd  = iim/2 -ndm )
55      PARAMETER( ip1jm  = iip1*jjm,  ip1jmp1= iip1*jjp1 )
56      PARAMETER( ip1jmi1= ip1jm - iip1 )
57      PARAMETER( ijp1llm= ip1jmp1 * llm, ijmllm= ip1jm * llm )
58      PARAMETER( mvar= ip1jmp1*( 2*llm+1) + ijmllm )
59      PARAMETER( jcfil=jjm/2+5, jcfllm=jcfil*llm )
60
61!-----------------------------------------------------------------------
62!
63! $Id: logic.h 1520 2011-05-23 11:37:09Z emillour $
64!
65!
66! NB: keep items of different kinds in seperate common blocs to avoid
67!     "misaligned commons" issues
68!-----------------------------------------------------------------------
69! INCLUDE 'logic.h'
70
71      COMMON/logicl/ purmats,forward,leapf,apphys,                      &
72     &  statcl,conser,apdiss,apdelq,saison,ecripar,fxyhypb,ysinus       &
73     &  ,read_start,ok_guide,ok_strato,tidal,ok_gradsfile               &
74     &  ,ok_limit,ok_etat0,hybrid                                       &
75     &  ,moyzon_mu,moyzon_ch
76
77      COMMON/logici/ iflag_phys,iflag_trac
78     
79      LOGICAL purmats,forward,leapf,apphys,statcl,conser,               &
80     & apdiss,apdelq,saison,ecripar,fxyhypb,ysinus                      &
81     &  ,read_start,ok_guide,ok_strato,tidal,ok_gradsfile               &
82     &  ,ok_limit,ok_etat0
83      logical hybrid ! vertical coordinate is hybrid if true (sigma otherwise)
84                     ! (only used if disvert_type==2)
85      logical moyzon_mu,moyzon_ch ! used for zonal averages in Titan
86
87      integer iflag_phys,iflag_trac
88!$OMP THREADPRIVATE(/logicl/)
89!$OMP THREADPRIVATE(/logici/)
90!-----------------------------------------------------------------------
91!
92! $Id: comvert.h 1654 2012-09-24 15:07:18Z aslmd $
93!
94!-----------------------------------------------------------------------
95!   INCLUDE 'comvert.h'
96
97      COMMON/comvertr/ap(llm+1),bp(llm+1),presnivs(llm),dpres(llm),     &
98     &               pa,preff,nivsigs(llm),nivsig(llm+1),               &
99     &               aps(llm),bps(llm),scaleheight,pseudoalt(llm)
100
101      common/comverti/disvert_type, pressure_exner
102
103      real ap     ! hybrid pressure contribution at interlayers
104      real bp     ! hybrid sigma contribution at interlayer
105      real presnivs ! (reference) pressure at mid-layers
106      real dpres
107      real pa     ! reference pressure (Pa) at which hybrid coordinates
108                  ! become purely pressure
109      real preff  ! reference surface pressure (Pa)
110      real nivsigs
111      real nivsig
112      real aps    ! hybrid pressure contribution at mid-layers
113      real bps    ! hybrid sigma contribution at mid-layers
114      real scaleheight ! atmospheric (reference) scale height (km)
115      real pseudoalt ! pseudo-altitude of model levels (km), based on presnivs(),
116                     ! preff and scaleheight
117
118      integer disvert_type ! type of vertical discretization:
119                           ! 1: Earth (default for planet_type==earth),
120                           !     automatic generation
121                           ! 2: Planets (default for planet_type!=earth),
122                           !     using 'z2sig.def' (or 'esasig.def) file
123
124      logical pressure_exner
125!     compute pressure inside layers using Exner function, else use mean
126!     of pressure values at interfaces
127
128 !-----------------------------------------------------------------------
129!
130! $Id: comconst.h 1437 2010-09-30 08:29:10Z emillour $
131!
132!-----------------------------------------------------------------------
133! INCLUDE comconst.h
134
135      COMMON/comconsti/im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl,          &
136     &                 iflag_top_bound,mode_top_bound
137      COMMON/comconstr/dtvr,daysec,                                     &
138     & pi,dtphys,dtdiss,rad,r,kappa,cotot,unsim,g,omeg                  &
139     & ,dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta        &
140     & ,dissip_pupstart  ,tau_top_bound,                                &
141     & daylen,molmass, ihf
142      COMMON/cpdetvenus/cpp,nu_venus,t0_venus
143
144      INTEGER im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl
145      REAL dtvr ! dynamical time step (in s)
146      REAL daysec !length (in s) of a standard day
147      REAL pi    ! something like 3.14159....
148      REAL dtphys ! (s) time step for the physics
149      REAL dtdiss ! (s) time step for the dissipation
150      REAL rad ! (m) radius of the planet
151      REAL r ! Reduced Gas constant r=R/mu
152             ! with R=8.31.. J.K-1.mol-1, mu: mol mass of atmosphere (kg/mol)
153      REAL cpp   ! Cp
154      REAL kappa ! kappa=R/Cp
155      REAL cotot
156      REAL unsim ! = 1./iim
157      REAL g ! (m/s2) gravity
158      REAL omeg ! (rad/s) rotation rate of the planet
159! Dissipation factors, for Earth model:
160      REAL dissip_factz,dissip_zref !dissip_deltaz
161! Dissipation factors, for other planets:
162      REAL dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta
163      REAL dissip_pupstart
164      INTEGER iflag_top_bound,mode_top_bound
165      REAL tau_top_bound
166      REAL daylen ! length of solar day, in 'standard' day length
167      REAL molmass ! (g/mol) molar mass of the atmosphere
168
169      REAL nu_venus,t0_venus ! coeffs needed for Cp(T), Venus atmosphere
170      REAL ihf  ! (W/m2) intrinsic heat flux for giant planets
171
172
173!-----------------------------------------------------------------------
174
175c
176c   Arguments:
177c   ----------
178      REAL masse(ip1jmp1,llm),pente_max
179c      REAL masse(iip1,jjp1,llm),pente_max
180      REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
181      REAL q(ip1jmp1,llm)
182c      REAL q(iip1,jjp1,llm)
183      REAL w(ip1jmp1,llm),pdt
184c
185c      Local 
186c   ---------
187c
188      INTEGER i,ij,l,j,ii
189      INTEGER ijlqmin,iqmin,jqmin,lqmin
190c
191      REAL zm(ip1jmp1,llm),newmasse
192      REAL mu(ip1jmp1,llm)
193      REAL mv(ip1jm,llm)
194      REAL mw(ip1jmp1,llm+1)
195      REAL zq(ip1jmp1,llm),zz
196      REAL dqx(ip1jmp1,llm),dqy(ip1jmp1,llm),dqz(ip1jmp1,llm)
197      REAL second,temps0,temps1,temps2,temps3
198      REAL ztemps1,ztemps2,ztemps3
199      REAL zzpbar, zzw
200      LOGICAL testcpu
201      SAVE testcpu
202      SAVE temps1,temps2,temps3
203      INTEGER iminn,imaxx
204
205      REAL qmin,qmax
206      DATA qmin,qmax/0.,1.e33/
207      DATA testcpu/.false./
208      DATA temps1,temps2,temps3/0.,0.,0./
209      INTEGER ijb,ije
210      type(request) :: MyRequest1
211      type(request) :: MyRequest2
212
213      call SetTag(MyRequest1,100)
214      call SetTag(MyRequest2,101)
215     
216      zzpbar = 0.5 * pdt
217      zzw    = pdt
218     
219      ijb=ij_begin
220      ije=ij_end
221      if (pole_nord) ijb=ijb+iip1
222      if (pole_sud)  ije=ije-iip1
223     
224      DO l=1,llm
225        DO ij = ijb,ije
226            mu(ij,l)=pbaru(ij,l) * zzpbar
227        ENDDO
228      ENDDO
229     
230      ijb=ij_begin-iip1
231      ije=ij_end
232      if (pole_nord) ijb=ij_begin
233      if (pole_sud)  ije=ij_end-iip1
234     
235      DO l=1,llm
236        DO ij=ijb,ije
237           mv(ij,l)=pbarv(ij,l) * zzpbar
238        ENDDO
239      ENDDO
240     
241      ijb=ij_begin
242      ije=ij_end
243     
244      DO l=1,llm
245        DO ij=ijb,ije
246           mw(ij,l)=w(ij,l) * zzw
247        ENDDO
248      ENDDO
249
250      DO ij=ijb,ije
251         mw(ij,llm+1)=0.
252      ENDDO
253     
254c      CALL SCOPY(ijp1llm,q,1,zq,1)
255c      CALL SCOPY(ijp1llm,masse,1,zm,1)
256       
257       ijb=ij_begin
258       ije=ij_end
259       zq(ijb:ije,:)=q(ijb:ije,:)
260       zm(ijb:ije,:)=masse(ijb:ije,:)
261     
262     
263c       print*,'Entree vlx1'
264c       call minmaxq(zq,qmin,qmax,'avant vlx     ')
265      call vlx_p(zq,pente_max,zm,mu,ij_begin,ij_begin+2*iip1-1)
266      call vlx_p(zq,pente_max,zm,mu,ij_end-2*iip1+1,ij_end)
267      call VTb(VTHallo)
268      call Register_Hallo(zq,ip1jmp1,llm,2,2,2,2,MyRequest1)
269      call Register_Hallo(zm,ip1jmp1,llm,1,1,1,1,MyRequest1)
270      call SendRequest(MyRequest1)
271      call VTe(VTHallo)
272      call vlx_p(zq,pente_max,zm,mu,ij_begin+2*iip1,ij_end-2*iip1)
273c      call vlx_p(zq,pente_max,zm,mu,ij_begin,ij_end)
274      call VTb(VTHallo)
275      call WaitRecvRequest(MyRequest1)
276      call VTe(VTHallo)
277
278     
279c       print*,'Sortie vlx1'
280c       call minmaxq(zq,qmin,qmax,'apres vlx1    ')
281
282c        print*,'Entree vly1'
283c      call exchange_hallo(zq,ip1jmp1,llm,2,2)
284c      call exchange_hallo(zm,ip1jmp1,llm,1,1)
285     
286      call vly_p(zq,pente_max,zm,mv)
287c       call minmaxq(zq,qmin,qmax,'apres vly1     ')
288c       print*,'Sortie vly1'
289      call vlz_p(zq,pente_max,zm,mw,ij_begin,ij_begin+2*iip1-1)
290      call vlz_p(zq,pente_max,zm,mw,ij_end-2*iip1+1,ij_end)
291      call VTb(VTHallo)
292      call Register_Hallo(zq,ip1jmp1,llm,2,2,2,2,MyRequest2)
293      call Register_Hallo(zm,ip1jmp1,llm,1,1,1,1,MyRequest2)
294      call SendRequest(MyRequest2)
295      call VTe(VTHallo)
296      call vlz_p(zq,pente_max,zm,mw,ij_begin+2*iip1,ij_end-2*iip1)
297      call VTb(VTHallo)
298      call WaitRecvRequest(MyRequest2)
299           
300      call VTe(VTHallo)
301     
302c       call minmaxq(zq,qmin,qmax,'apres vlz     ')
303
304
305     
306     
307      call vly_p(zq,pente_max,zm,mv)
308c       call minmaxq(zq,qmin,qmax,'apres vly     ')
309
310
311      call vlx_p(zq,pente_max,zm,mu,ij_begin,ij_end)
312c       call minmaxq(zq,qmin,qmax,'apres vlx2    ')
313
314       
315      ijb=ij_begin
316      ije=ij_end
317       
318      DO l=1,llm
319         DO ij=ijb,ije
320           q(ij,l)=zq(ij,l)
321         ENDDO
322      ENDDO
323     
324     
325      DO l=1,llm
326         DO ij=ijb,ije-iip1+1,iip1
327            q(ij+iim,l)=q(ij,l)
328         ENDDO
329      ENDDO
330
331      call WaitSendRequest(MyRequest1) 
332      call WaitSendRequest(MyRequest2)
333     
334      RETURN
335      END
336     
337     
338      SUBROUTINE vlx_p(q,pente_max,masse,u_m,ijb_x,ije_x)
339
340c     Auteurs:   P.Le Van, F.Hourdin, F.Forget 
341c
342c    ********************************************************************
343c     Shema  d'advection " pseudo amont " .
344c    ********************************************************************
345c     nq,iq,q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....
346c
347c
348c   --------------------------------------------------------------------
349      USE Parallel_lmdz
350      IMPLICIT NONE
351c
352!-----------------------------------------------------------------------
353!   INCLUDE 'dimensions.h'
354!
355!   dimensions.h contient les dimensions du modele
356!   ndm est tel que iim=2**ndm
357!-----------------------------------------------------------------------
358
359      INTEGER iim,jjm,llm,ndm
360
361      PARAMETER (iim= 128,jjm=96,llm=64,ndm=1)
362
363!-----------------------------------------------------------------------
364!
365! $Header$
366!
367!
368!  ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
369!                 veillez  n'utiliser que des ! pour les commentaires
370!                 et  bien positionner les & des lignes de continuation
371!                 (les placer en colonne 6 et en colonne 73)
372!
373!
374!-----------------------------------------------------------------------
375!   INCLUDE 'paramet.h'
376
377      INTEGER  iip1,iip2,iip3,jjp1,llmp1,llmp2,llmm1
378      INTEGER  kftd,ip1jm,ip1jmp1,ip1jmi1,ijp1llm
379      INTEGER  ijmllm,mvar
380      INTEGER jcfil,jcfllm
381
382      PARAMETER( iip1= iim+1,iip2=iim+2,iip3=iim+3                       &
383     &    ,jjp1=jjm+1-1/jjm)
384      PARAMETER( llmp1 = llm+1,  llmp2 = llm+2, llmm1 = llm-1 )
385      PARAMETER( kftd  = iim/2 -ndm )
386      PARAMETER( ip1jm  = iip1*jjm,  ip1jmp1= iip1*jjp1 )
387      PARAMETER( ip1jmi1= ip1jm - iip1 )
388      PARAMETER( ijp1llm= ip1jmp1 * llm, ijmllm= ip1jm * llm )
389      PARAMETER( mvar= ip1jmp1*( 2*llm+1) + ijmllm )
390      PARAMETER( jcfil=jjm/2+5, jcfllm=jcfil*llm )
391
392!-----------------------------------------------------------------------
393!
394! $Id: logic.h 1520 2011-05-23 11:37:09Z emillour $
395!
396!
397! NB: keep items of different kinds in seperate common blocs to avoid
398!     "misaligned commons" issues
399!-----------------------------------------------------------------------
400! INCLUDE 'logic.h'
401
402      COMMON/logicl/ purmats,forward,leapf,apphys,                      &
403     &  statcl,conser,apdiss,apdelq,saison,ecripar,fxyhypb,ysinus       &
404     &  ,read_start,ok_guide,ok_strato,tidal,ok_gradsfile               &
405     &  ,ok_limit,ok_etat0,hybrid                                       &
406     &  ,moyzon_mu,moyzon_ch
407
408      COMMON/logici/ iflag_phys,iflag_trac
409     
410      LOGICAL purmats,forward,leapf,apphys,statcl,conser,               &
411     & apdiss,apdelq,saison,ecripar,fxyhypb,ysinus                      &
412     &  ,read_start,ok_guide,ok_strato,tidal,ok_gradsfile               &
413     &  ,ok_limit,ok_etat0
414      logical hybrid ! vertical coordinate is hybrid if true (sigma otherwise)
415                     ! (only used if disvert_type==2)
416      logical moyzon_mu,moyzon_ch ! used for zonal averages in Titan
417
418      integer iflag_phys,iflag_trac
419!$OMP THREADPRIVATE(/logicl/)
420!$OMP THREADPRIVATE(/logici/)
421!-----------------------------------------------------------------------
422!
423! $Id: comvert.h 1654 2012-09-24 15:07:18Z aslmd $
424!
425!-----------------------------------------------------------------------
426!   INCLUDE 'comvert.h'
427
428      COMMON/comvertr/ap(llm+1),bp(llm+1),presnivs(llm),dpres(llm),     &
429     &               pa,preff,nivsigs(llm),nivsig(llm+1),               &
430     &               aps(llm),bps(llm),scaleheight,pseudoalt(llm)
431
432      common/comverti/disvert_type, pressure_exner
433
434      real ap     ! hybrid pressure contribution at interlayers
435      real bp     ! hybrid sigma contribution at interlayer
436      real presnivs ! (reference) pressure at mid-layers
437      real dpres
438      real pa     ! reference pressure (Pa) at which hybrid coordinates
439                  ! become purely pressure
440      real preff  ! reference surface pressure (Pa)
441      real nivsigs
442      real nivsig
443      real aps    ! hybrid pressure contribution at mid-layers
444      real bps    ! hybrid sigma contribution at mid-layers
445      real scaleheight ! atmospheric (reference) scale height (km)
446      real pseudoalt ! pseudo-altitude of model levels (km), based on presnivs(),
447                     ! preff and scaleheight
448
449      integer disvert_type ! type of vertical discretization:
450                           ! 1: Earth (default for planet_type==earth),
451                           !     automatic generation
452                           ! 2: Planets (default for planet_type!=earth),
453                           !     using 'z2sig.def' (or 'esasig.def) file
454
455      logical pressure_exner
456!     compute pressure inside layers using Exner function, else use mean
457!     of pressure values at interfaces
458
459 !-----------------------------------------------------------------------
460!
461! $Id: comconst.h 1437 2010-09-30 08:29:10Z emillour $
462!
463!-----------------------------------------------------------------------
464! INCLUDE comconst.h
465
466      COMMON/comconsti/im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl,          &
467     &                 iflag_top_bound,mode_top_bound
468      COMMON/comconstr/dtvr,daysec,                                     &
469     & pi,dtphys,dtdiss,rad,r,kappa,cotot,unsim,g,omeg                  &
470     & ,dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta        &
471     & ,dissip_pupstart  ,tau_top_bound,                                &
472     & daylen,molmass, ihf
473      COMMON/cpdetvenus/cpp,nu_venus,t0_venus
474
475      INTEGER im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl
476      REAL dtvr ! dynamical time step (in s)
477      REAL daysec !length (in s) of a standard day
478      REAL pi    ! something like 3.14159....
479      REAL dtphys ! (s) time step for the physics
480      REAL dtdiss ! (s) time step for the dissipation
481      REAL rad ! (m) radius of the planet
482      REAL r ! Reduced Gas constant r=R/mu
483             ! with R=8.31.. J.K-1.mol-1, mu: mol mass of atmosphere (kg/mol)
484      REAL cpp   ! Cp
485      REAL kappa ! kappa=R/Cp
486      REAL cotot
487      REAL unsim ! = 1./iim
488      REAL g ! (m/s2) gravity
489      REAL omeg ! (rad/s) rotation rate of the planet
490! Dissipation factors, for Earth model:
491      REAL dissip_factz,dissip_zref !dissip_deltaz
492! Dissipation factors, for other planets:
493      REAL dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta
494      REAL dissip_pupstart
495      INTEGER iflag_top_bound,mode_top_bound
496      REAL tau_top_bound
497      REAL daylen ! length of solar day, in 'standard' day length
498      REAL molmass ! (g/mol) molar mass of the atmosphere
499
500      REAL nu_venus,t0_venus ! coeffs needed for Cp(T), Venus atmosphere
501      REAL ihf  ! (W/m2) intrinsic heat flux for giant planets
502
503
504!-----------------------------------------------------------------------
505c
506c
507c   Arguments:
508c   ----------
509      REAL masse(ip1jmp1,llm),pente_max
510      REAL u_m( ip1jmp1,llm ),pbarv( iip1,jjm,llm)
511      REAL q(ip1jmp1,llm)
512      REAL w(ip1jmp1,llm)
513c
514c      Local 
515c   ---------
516c
517      INTEGER ij,l,j,i,iju,ijq,indu(ip1jmp1),niju
518      INTEGER n0,iadvplus(ip1jmp1,llm),nl(llm)
519c
520      REAL new_m,zu_m,zdum(ip1jmp1,llm)
521      REAL sigu(ip1jmp1),dxq(ip1jmp1,llm),dxqu(ip1jmp1)
522      REAL zz(ip1jmp1)
523      REAL adxqu(ip1jmp1),dxqmax(ip1jmp1,llm)
524      REAL u_mq(ip1jmp1,llm)
525
526      Logical extremum
527
528      REAL      SSUM
529      EXTERNAL  SSUM
530
531      REAL z1,z2,z3
532
533      INTEGER ijb,ije,ijb_x,ije_x
534     
535c   calcul de la pente a droite et a gauche de la maille
536
537      ijb=ijb_x
538      ije=ije_x
539       
540      if (pole_nord.and.ijb==1) ijb=ijb+iip1
541      if (pole_sud.and.ije==ip1jmp1)  ije=ije-iip1
542         
543      IF (pente_max.gt.-1.e-5) THEN
544c       IF (pente_max.gt.10) THEN
545
546c   calcul des pentes avec limitation, Van Leer scheme I:
547c   -----------------------------------------------------
548
549c   calcul de la pente aux points u
550c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)         
551         DO l = 1, llm
552           
553            DO ij=ijb,ije-1
554               dxqu(ij)=q(ij+1,l)-q(ij,l)
555c              IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
556c              sigu(ij)=u_m(ij,l)/masse(ij,l)
557            ENDDO
558            DO ij=ijb+iip1-1,ije,iip1
559               dxqu(ij)=dxqu(ij-iim)
560c              sigu(ij)=sigu(ij-iim)
561            ENDDO
562
563            DO ij=ijb,ije
564               adxqu(ij)=abs(dxqu(ij))
565            ENDDO
566
567c   calcul de la pente maximum dans la maille en valeur absolue
568
569            DO ij=ijb+1,ije
570               dxqmax(ij,l)=pente_max*
571     ,      min(adxqu(ij-1),adxqu(ij))
572c limitation subtile
573c    ,      min(adxqu(ij-1)/sigu(ij-1),adxqu(ij)/(1.-sigu(ij)))
574         
575
576            ENDDO
577
578            DO ij=ijb+iip1-1,ije,iip1
579               dxqmax(ij-iim,l)=dxqmax(ij,l)
580            ENDDO
581
582            DO ij=ijb+1,ije
583               IF(dxqu(ij-1)*dxqu(ij).gt.0) THEN
584                  dxq(ij,l)=dxqu(ij-1)+dxqu(ij)
585               ELSE
586c   extremum local
587                  dxq(ij,l)=0.
588               ENDIF
589               dxq(ij,l)=0.5*dxq(ij,l)
590               dxq(ij,l)=
591     ,         sign(min(abs(dxq(ij,l)),dxqmax(ij,l)),dxq(ij,l))
592            ENDDO
593
594         ENDDO ! l=1,llm
595c$OMP END DO NOWAIT
596c       print*,'Ok calcul des pentes'
597
598      ELSE ! (pente_max.lt.-1.e-5)
599
600c   Pentes produits:
601c   ----------------
602c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
603         DO l = 1, llm
604            DO ij=ijb,ije-1
605               dxqu(ij)=q(ij+1,l)-q(ij,l)
606            ENDDO
607            DO ij=ijb+iip1-1,ije,iip1
608               dxqu(ij)=dxqu(ij-iim)
609            ENDDO
610
611            DO ij=ijb+1,ije
612               zz(ij)=dxqu(ij-1)*dxqu(ij)
613               zz(ij)=zz(ij)+zz(ij)
614               IF(zz(ij).gt.0) THEN
615                  dxq(ij,l)=zz(ij)/(dxqu(ij-1)+dxqu(ij))
616               ELSE
617c   extremum local
618                  dxq(ij,l)=0.
619               ENDIF
620            ENDDO
621
622         ENDDO
623c$OMP END DO NOWAIT
624      ENDIF ! (pente_max.lt.-1.e-5)
625
626c   bouclage de la pente en iip1:
627c   -----------------------------
628c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
629      DO l=1,llm
630         DO ij=ijb+iip1-1,ije,iip1
631            dxq(ij-iim,l)=dxq(ij,l)
632         ENDDO
633         DO ij=ijb,ije
634            iadvplus(ij,l)=0
635         ENDDO
636
637      ENDDO
638c$OMP END DO NOWAIT
639c        print*,'Bouclage en iip1'
640
641c   calcul des flux a gauche et a droite
642
643c   on cumule le flux correspondant a toutes les mailles dont la masse
644c   au travers de la paroi pENDant le pas de temps.
645c       print*,'Cumule ....'
646c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
647      DO l=1,llm
648       DO ij=ijb,ije-1
649c       print*,'masse(',ij,')=',masse(ij,l)
650          IF (u_m(ij,l).gt.0.) THEN
651             zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l)
652             u_mq(ij,l)=u_m(ij,l)*(q(ij,l)+0.5*zdum(ij,l)*dxq(ij,l))
653          ELSE
654             zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l)
655             u_mq(ij,l)=u_m(ij,l)*(q(ij+1,l)-0.5*zdum(ij,l)*dxq(ij+1,l))
656          ENDIF
657       ENDDO
658      ENDDO
659c$OMP END DO NOWAIT
660c       stop
661
662c       go to 9999
663c   detection des points ou on advecte plus que la masse de la
664c   maille
665c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
666      DO l=1,llm
667         DO ij=ijb,ije-1
668            IF(zdum(ij,l).lt.0) THEN
669               iadvplus(ij,l)=1
670               u_mq(ij,l)=0.
671            ENDIF
672         ENDDO
673      ENDDO
674c$OMP END DO NOWAIT
675c       print*,'Ok test 1'
676c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
677      DO l=1,llm
678       DO ij=ijb+iip1-1,ije,iip1
679          iadvplus(ij,l)=iadvplus(ij-iim,l)
680       ENDDO
681      ENDDO
682c$OMP END DO NOWAIT
683c        print*,'Ok test 2'
684
685
686c   traitement special pour le cas ou on advecte en longitude plus que le
687c   contenu de la maille.
688c   cette partie est mal vectorisee.
689
690c  calcul du nombre de maille sur lequel on advecte plus que la maille.
691
692      n0=0
693c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
694      DO l=1,llm
695         nl(l)=0
696         DO ij=ijb,ije
697            nl(l)=nl(l)+iadvplus(ij,l)
698         ENDDO
699         n0=n0+nl(l)
700      ENDDO
701c$OMP END DO NOWAIT
702cym      IF(n0.gt.1) THEN
703cym      IF(n0.gt.0) THEN
704
705c      PRINT*,'Nombre de points pour lesquels on advect plus que le'
706c     &       ,'contenu de la maille : ',n0
707c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
708         DO l=1,llm
709            IF(nl(l).gt.0) THEN
710               iju=0
711c   indicage des mailles concernees par le traitement special
712               DO ij=ijb,ije
713                  IF(iadvplus(ij,l).eq.1.and.mod(ij,iip1).ne.0) THEN
714                     iju=iju+1
715                     indu(iju)=ij
716                  ENDIF
717               ENDDO
718               niju=iju
719c              PRINT*,'niju,nl',niju,nl(l)
720
721c  traitement des mailles
722               DO iju=1,niju
723                  ij=indu(iju)
724                  j=(ij-1)/iip1+1
725                  zu_m=u_m(ij,l)
726                  u_mq(ij,l)=0.
727                  IF(zu_m.gt.0.) THEN
728                     ijq=ij
729                     i=ijq-(j-1)*iip1
730c   accumulation pour les mailles completements advectees
731                     do while(zu_m.gt.masse(ijq,l))
732                        u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l)
733                        zu_m=zu_m-masse(ijq,l)
734                        i=mod(i-2+iim,iim)+1
735                        ijq=(j-1)*iip1+i
736                     ENDDO
737c   ajout de la maille non completement advectee
738                     u_mq(ij,l)=u_mq(ij,l)+zu_m*
739     &               (q(ijq,l)+0.5*(1.-zu_m/masse(ijq,l))*dxq(ijq,l))
740                  ELSE
741                     ijq=ij+1
742                     i=ijq-(j-1)*iip1
743c   accumulation pour les mailles completements advectees
744                     do while(-zu_m.gt.masse(ijq,l))
745                        u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l)
746                        zu_m=zu_m+masse(ijq,l)
747                        i=mod(i,iim)+1
748                        ijq=(j-1)*iip1+i
749                     ENDDO
750c   ajout de la maille non completement advectee
751                     u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l)-
752     &               0.5*(1.+zu_m/masse(ijq,l))*dxq(ijq,l))
753                  ENDIF
754               ENDDO
755            ENDIF
756         ENDDO
757c$OMP END DO NOWAIT
758cym      ENDIF  ! n0.gt.0
7599999    continue
760
761
762c   bouclage en latitude
763c       print*,'Avant bouclage en latitude'
764c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
765      DO l=1,llm
766        DO ij=ijb+iip1-1,ije,iip1
767           u_mq(ij,l)=u_mq(ij-iim,l)
768        ENDDO
769      ENDDO
770c$OMP END DO NOWAIT
771
772c   calcul des tENDances
773c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
774      DO l=1,llm
775         DO ij=ijb+1,ije
776            new_m=masse(ij,l)+u_m(ij-1,l)-u_m(ij,l)
777            q(ij,l)=(q(ij,l)*masse(ij,l)+
778     &      u_mq(ij-1,l)-u_mq(ij,l))
779     &      /new_m
780            masse(ij,l)=new_m
781         ENDDO
782c   ModIF Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
783         DO ij=ijb+iip1-1,ije,iip1
784            q(ij-iim,l)=q(ij,l)
785            masse(ij-iim,l)=masse(ij,l)
786         ENDDO
787      ENDDO
788c$OMP END DO NOWAIT
789c     CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
790c     CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1)
791
792
793      RETURN
794      END
795
796
797      SUBROUTINE vly_p(q,pente_max,masse,masse_adv_v)
798c
799c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
800c
801c    ********************************************************************
802c     Shema  d'advection " pseudo amont " .
803c    ********************************************************************
804c     q,masse_adv_v,w sont des arguments d'entree  pour le s-pg ....
805c     dq               sont des arguments de sortie pour le s-pg ....
806c
807c
808c   --------------------------------------------------------------------
809      USE parallel_lmdz
810      IMPLICIT NONE
811c
812!-----------------------------------------------------------------------
813!   INCLUDE 'dimensions.h'
814!
815!   dimensions.h contient les dimensions du modele
816!   ndm est tel que iim=2**ndm
817!-----------------------------------------------------------------------
818
819      INTEGER iim,jjm,llm,ndm
820
821      PARAMETER (iim= 128,jjm=96,llm=64,ndm=1)
822
823!-----------------------------------------------------------------------
824!
825! $Header$
826!
827!
828!  ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
829!                 veillez  n'utiliser que des ! pour les commentaires
830!                 et  bien positionner les & des lignes de continuation
831!                 (les placer en colonne 6 et en colonne 73)
832!
833!
834!-----------------------------------------------------------------------
835!   INCLUDE 'paramet.h'
836
837      INTEGER  iip1,iip2,iip3,jjp1,llmp1,llmp2,llmm1
838      INTEGER  kftd,ip1jm,ip1jmp1,ip1jmi1,ijp1llm
839      INTEGER  ijmllm,mvar
840      INTEGER jcfil,jcfllm
841
842      PARAMETER( iip1= iim+1,iip2=iim+2,iip3=iim+3                       &
843     &    ,jjp1=jjm+1-1/jjm)
844      PARAMETER( llmp1 = llm+1,  llmp2 = llm+2, llmm1 = llm-1 )
845      PARAMETER( kftd  = iim/2 -ndm )
846      PARAMETER( ip1jm  = iip1*jjm,  ip1jmp1= iip1*jjp1 )
847      PARAMETER( ip1jmi1= ip1jm - iip1 )
848      PARAMETER( ijp1llm= ip1jmp1 * llm, ijmllm= ip1jm * llm )
849      PARAMETER( mvar= ip1jmp1*( 2*llm+1) + ijmllm )
850      PARAMETER( jcfil=jjm/2+5, jcfllm=jcfil*llm )
851
852!-----------------------------------------------------------------------
853!
854! $Id: logic.h 1520 2011-05-23 11:37:09Z emillour $
855!
856!
857! NB: keep items of different kinds in seperate common blocs to avoid
858!     "misaligned commons" issues
859!-----------------------------------------------------------------------
860! INCLUDE 'logic.h'
861
862      COMMON/logicl/ purmats,forward,leapf,apphys,                      &
863     &  statcl,conser,apdiss,apdelq,saison,ecripar,fxyhypb,ysinus       &
864     &  ,read_start,ok_guide,ok_strato,tidal,ok_gradsfile               &
865     &  ,ok_limit,ok_etat0,hybrid                                       &
866     &  ,moyzon_mu,moyzon_ch
867
868      COMMON/logici/ iflag_phys,iflag_trac
869     
870      LOGICAL purmats,forward,leapf,apphys,statcl,conser,               &
871     & apdiss,apdelq,saison,ecripar,fxyhypb,ysinus                      &
872     &  ,read_start,ok_guide,ok_strato,tidal,ok_gradsfile               &
873     &  ,ok_limit,ok_etat0
874      logical hybrid ! vertical coordinate is hybrid if true (sigma otherwise)
875                     ! (only used if disvert_type==2)
876      logical moyzon_mu,moyzon_ch ! used for zonal averages in Titan
877
878      integer iflag_phys,iflag_trac
879!$OMP THREADPRIVATE(/logicl/)
880!$OMP THREADPRIVATE(/logici/)
881!-----------------------------------------------------------------------
882!
883! $Id: comvert.h 1654 2012-09-24 15:07:18Z aslmd $
884!
885!-----------------------------------------------------------------------
886!   INCLUDE 'comvert.h'
887
888      COMMON/comvertr/ap(llm+1),bp(llm+1),presnivs(llm),dpres(llm),     &
889     &               pa,preff,nivsigs(llm),nivsig(llm+1),               &
890     &               aps(llm),bps(llm),scaleheight,pseudoalt(llm)
891
892      common/comverti/disvert_type, pressure_exner
893
894      real ap     ! hybrid pressure contribution at interlayers
895      real bp     ! hybrid sigma contribution at interlayer
896      real presnivs ! (reference) pressure at mid-layers
897      real dpres
898      real pa     ! reference pressure (Pa) at which hybrid coordinates
899                  ! become purely pressure
900      real preff  ! reference surface pressure (Pa)
901      real nivsigs
902      real nivsig
903      real aps    ! hybrid pressure contribution at mid-layers
904      real bps    ! hybrid sigma contribution at mid-layers
905      real scaleheight ! atmospheric (reference) scale height (km)
906      real pseudoalt ! pseudo-altitude of model levels (km), based on presnivs(),
907                     ! preff and scaleheight
908
909      integer disvert_type ! type of vertical discretization:
910                           ! 1: Earth (default for planet_type==earth),
911                           !     automatic generation
912                           ! 2: Planets (default for planet_type!=earth),
913                           !     using 'z2sig.def' (or 'esasig.def) file
914
915      logical pressure_exner
916!     compute pressure inside layers using Exner function, else use mean
917!     of pressure values at interfaces
918
919 !-----------------------------------------------------------------------
920!
921! $Id: comconst.h 1437 2010-09-30 08:29:10Z emillour $
922!
923!-----------------------------------------------------------------------
924! INCLUDE comconst.h
925
926      COMMON/comconsti/im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl,          &
927     &                 iflag_top_bound,mode_top_bound
928      COMMON/comconstr/dtvr,daysec,                                     &
929     & pi,dtphys,dtdiss,rad,r,kappa,cotot,unsim,g,omeg                  &
930     & ,dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta        &
931     & ,dissip_pupstart  ,tau_top_bound,                                &
932     & daylen,molmass, ihf
933      COMMON/cpdetvenus/cpp,nu_venus,t0_venus
934
935      INTEGER im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl
936      REAL dtvr ! dynamical time step (in s)
937      REAL daysec !length (in s) of a standard day
938      REAL pi    ! something like 3.14159....
939      REAL dtphys ! (s) time step for the physics
940      REAL dtdiss ! (s) time step for the dissipation
941      REAL rad ! (m) radius of the planet
942      REAL r ! Reduced Gas constant r=R/mu
943             ! with R=8.31.. J.K-1.mol-1, mu: mol mass of atmosphere (kg/mol)
944      REAL cpp   ! Cp
945      REAL kappa ! kappa=R/Cp
946      REAL cotot
947      REAL unsim ! = 1./iim
948      REAL g ! (m/s2) gravity
949      REAL omeg ! (rad/s) rotation rate of the planet
950! Dissipation factors, for Earth model:
951      REAL dissip_factz,dissip_zref !dissip_deltaz
952! Dissipation factors, for other planets:
953      REAL dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta
954      REAL dissip_pupstart
955      INTEGER iflag_top_bound,mode_top_bound
956      REAL tau_top_bound
957      REAL daylen ! length of solar day, in 'standard' day length
958      REAL molmass ! (g/mol) molar mass of the atmosphere
959
960      REAL nu_venus,t0_venus ! coeffs needed for Cp(T), Venus atmosphere
961      REAL ihf  ! (W/m2) intrinsic heat flux for giant planets
962
963
964!-----------------------------------------------------------------------
965!
966! $Header$
967!
968!CDK comgeom
969      COMMON/comgeom/                                                   &
970     & cu(ip1jmp1),cv(ip1jm),unscu2(ip1jmp1),unscv2(ip1jm),             &
971     & aire(ip1jmp1),airesurg(ip1jmp1),aireu(ip1jmp1),                  &
972     & airev(ip1jm),unsaire(ip1jmp1),apoln,apols,                       &
973     & unsairez(ip1jm),airuscv2(ip1jm),airvscu2(ip1jm),                 &
974     & aireij1(ip1jmp1),aireij2(ip1jmp1),aireij3(ip1jmp1),              &
975     & aireij4(ip1jmp1),alpha1(ip1jmp1),alpha2(ip1jmp1),                &
976     & alpha3(ip1jmp1),alpha4(ip1jmp1),alpha1p2(ip1jmp1),               &
977     & alpha1p4(ip1jmp1),alpha2p3(ip1jmp1),alpha3p4(ip1jmp1),           &
978     & fext(ip1jm),constang(ip1jmp1),rlatu(jjp1),rlatv(jjm),            &
979     & rlonu(iip1),rlonv(iip1),cuvsurcv(ip1jm),cvsurcuv(ip1jm),         &
980     & cvusurcu(ip1jmp1),cusurcvu(ip1jmp1),cuvscvgam1(ip1jm),           &
981     & cuvscvgam2(ip1jm),cvuscugam1(ip1jmp1),                           &
982     & cvuscugam2(ip1jmp1),cvscuvgam(ip1jm),cuscvugam(ip1jmp1),         &
983     & unsapolnga1,unsapolnga2,unsapolsga1,unsapolsga2,                 &
984     & unsair_gam1(ip1jmp1),unsair_gam2(ip1jmp1),unsairz_gam(ip1jm),    &
985     & aivscu2gam(ip1jm),aiuscv2gam(ip1jm),xprimu(iip1),xprimv(iip1)
986
987!
988        REAL                                                            &
989     & cu,cv,unscu2,unscv2,aire,airesurg,aireu,airev,unsaire,apoln     ,&
990     & apols,unsairez,airuscv2,airvscu2,aireij1,aireij2,aireij3,aireij4,&
991     & alpha1,alpha2,alpha3,alpha4,alpha1p2,alpha1p4,alpha2p3,alpha3p4 ,&
992     & fext,constang,rlatu,rlatv,rlonu,rlonv,cuvscvgam1,cuvscvgam2     ,&
993     & cvuscugam1,cvuscugam2,cvscuvgam,cuscvugam,unsapolnga1,unsapolnga2&
994     & ,unsapolsga1,unsapolsga2,unsair_gam1,unsair_gam2,unsairz_gam    ,&
995     & aivscu2gam ,aiuscv2gam,cuvsurcv,cvsurcuv,cvusurcu,cusurcvu,xprimu&
996     & , xprimv
997!
998c
999c
1000c   Arguments:
1001c   ----------
1002      REAL masse(ip1jmp1,llm),pente_max
1003      REAL masse_adv_v( ip1jm,llm)
1004      REAL q(ip1jmp1,llm), dq( ip1jmp1,llm)
1005c
1006c      Local 
1007c   ---------
1008c
1009      INTEGER i,ij,l
1010c
1011      REAL airej2,airejjm,airescb(iim),airesch(iim)
1012      REAL dyq(ip1jmp1,llm),dyqv(ip1jm),zdvm(ip1jmp1,llm)
1013      REAL adyqv(ip1jm),dyqmax(ip1jmp1)
1014      REAL qbyv(ip1jm,llm)
1015
1016      REAL qpns,qpsn,appn,apps,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
1017c     REAL newq,oldmasse
1018      Logical extremum,first,testcpu
1019      REAL temps0,temps1,temps2,temps3,temps4,temps5,second
1020      SAVE temps0,temps1,temps2,temps3,temps4,temps5
1021c$OMP THREADPRIVATE(temps0,temps1,temps2,temps3,temps4,temps5)
1022      SAVE first,testcpu
1023c$OMP THREADPRIVATE(first,testcpu)
1024
1025      REAL convpn,convps,convmpn,convmps
1026      real massepn,masseps,qpn,qps
1027      REAL sinlon(iip1),sinlondlon(iip1)
1028      REAL coslon(iip1),coslondlon(iip1)
1029      SAVE sinlon,coslon,sinlondlon,coslondlon
1030c$OMP THREADPRIVATE(sinlon,coslon,sinlondlon,coslondlon)
1031      SAVE airej2,airejjm
1032c$OMP THREADPRIVATE(airej2,airejjm)
1033c
1034c
1035      REAL      SSUM
1036      EXTERNAL  SSUM
1037
1038      DATA first,testcpu/.true.,.false./
1039      DATA temps0,temps1,temps2,temps3,temps4,temps5/0.,0.,0.,0.,0.,0./
1040      INTEGER ijb,ije
1041
1042      IF(first) THEN
1043c         PRINT*,'Shema  Amont nouveau  appele dans  Vanleer   '
1044         first=.false.
1045         do i=2,iip1
1046            coslon(i)=cos(rlonv(i))
1047            sinlon(i)=sin(rlonv(i))
1048            coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi
1049            sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
1050         ENDDO
1051         coslon(1)=coslon(iip1)
1052         coslondlon(1)=coslondlon(iip1)
1053         sinlon(1)=sinlon(iip1)
1054         sinlondlon(1)=sinlondlon(iip1)
1055         airej2 = SSUM( iim, aire(iip2), 1 )
1056         airejjm= SSUM( iim, aire(ip1jm -iim), 1 ) 
1057      ENDIF
1058
1059c
1060c       PRINT*,'CALCUL EN LATITUDE'
1061
1062c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
1063      DO l = 1, llm
1064c
1065c   --------------------------------
1066c      CALCUL EN LATITUDE
1067c   --------------------------------
1068
1069c   On commence par calculer la valeur du traceur moyenne sur le premier cercle
1070c   de latitude autour du pole (qpns pour le pole nord et qpsn pour
1071c    le pole nord) qui sera utilisee pour evaluer les pentes au pole.
1072     
1073      if (pole_nord) then
1074        DO i = 1, iim
1075          airescb(i) = aire(i+ iip1) * q(i+ iip1,l)
1076        ENDDO
1077        qpns   = SSUM( iim,  airescb ,1 ) / airej2
1078      endif
1079     
1080      if (pole_sud) then
1081        DO i = 1, iim
1082          airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l)
1083        ENDDO
1084        qpsn   = SSUM( iim,  airesch ,1 ) / airejjm
1085      endif
1086     
1087     
1088
1089c   calcul des pentes aux points v
1090
1091      ijb=ij_begin-2*iip1
1092      ije=ij_end+iip1
1093      if (pole_nord) ijb=ij_begin
1094      if (pole_sud)  ije=ij_end-iip1
1095     
1096      DO ij=ijb,ije
1097         dyqv(ij)=q(ij,l)-q(ij+iip1,l)
1098         adyqv(ij)=abs(dyqv(ij))
1099      ENDDO
1100
1101c   calcul des pentes aux points scalaires
1102      ijb=ij_begin-iip1
1103      ije=ij_end+iip1
1104      if (pole_nord) ijb=ij_begin+iip1
1105      if (pole_sud)  ije=ij_end-iip1
1106     
1107      DO ij=ijb,ije
1108         dyq(ij,l)=.5*(dyqv(ij-iip1)+dyqv(ij))
1109         dyqmax(ij)=min(adyqv(ij-iip1),adyqv(ij))
1110         dyqmax(ij)=pente_max*dyqmax(ij)
1111      ENDDO
1112
1113c   calcul des pentes aux poles
1114      IF (pole_nord) THEN
1115        DO ij=1,iip1
1116           dyq(ij,l)=qpns-q(ij+iip1,l)
1117        ENDDO
1118       
1119        dyn1=0.
1120        dyn2=0.
1121        DO ij=1,iim
1122          dyn1=dyn1+sinlondlon(ij)*dyq(ij,l)
1123          dyn2=dyn2+coslondlon(ij)*dyq(ij,l)
1124        ENDDO
1125        DO ij=1,iip1
1126          dyq(ij,l)=dyn1*sinlon(ij)+dyn2*coslon(ij)
1127        ENDDO
1128       
1129        DO ij=1,iip1
1130         dyq(ij,l)=0.
1131        ENDDO
1132c ym tout cela ne sert pas a grand chose
1133      ENDIF
1134     
1135      IF (pole_sud) THEN
1136
1137        DO ij=1,iip1
1138           dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l)-qpsn
1139        ENDDO
1140
1141        dys1=0.
1142        dys2=0.
1143
1144        DO ij=1,iim
1145          dys1=dys1+sinlondlon(ij)*dyq(ip1jm+ij,l)
1146          dys2=dys2+coslondlon(ij)*dyq(ip1jm+ij,l)
1147        ENDDO
1148
1149        DO ij=1,iip1
1150          dyq(ip1jm+ij,l)=dys1*sinlon(ij)+dys2*coslon(ij)
1151        ENDDO
1152       
1153        DO ij=1,iip1
1154         dyq(ip1jm+ij,l)=0.
1155        ENDDO
1156c ym tout cela ne sert pas a grand chose
1157      ENDIF
1158
1159c   filtrage de la derivee
1160
1161c   calcul des pentes limites aux poles
1162c ym partie inutile
1163c      goto 8888
1164c      fn=1.
1165c      fs=1.
1166c      DO ij=1,iim
1167c         IF(pente_max*adyqv(ij).lt.abs(dyq(ij,l))) THEN
1168c            fn=min(pente_max*adyqv(ij)/abs(dyq(ij,l)),fn)
1169c         ENDIF
1170c      IF(pente_max*adyqv(ij+ip1jm-iip1).lt.abs(dyq(ij+ip1jm,l))) THEN
1171c         fs=min(pente_max*adyqv(ij+ip1jm-iip1)/abs(dyq(ij+ip1jm,l)),fs)
1172c         ENDIF
1173c      ENDDO
1174c      DO ij=1,iip1
1175c         dyq(ij,l)=fn*dyq(ij,l)
1176c         dyq(ip1jm+ij,l)=fs*dyq(ip1jm+ij,l)
1177c      ENDDO
1178c 8888    continue
1179
1180
1181CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1182C  En memoire de dIFferents tests sur la 
1183C  limitation des pentes aux poles.
1184CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1185C     PRINT*,dyq(1)
1186C     PRINT*,dyqv(iip1+1)
1187C     appn=abs(dyq(1)/dyqv(iip1+1))
1188C     PRINT*,dyq(ip1jm+1)
1189C     PRINT*,dyqv(ip1jm-iip1+1)
1190C     apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
1191C     DO ij=2,iim
1192C        appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
1193C        apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
1194C     ENDDO
1195C     appn=min(pente_max/appn,1.)
1196C     apps=min(pente_max/apps,1.)
1197C
1198C
1199C   cas ou on a un extremum au pole
1200C
1201C     IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
1202C    &   appn=0.
1203C     IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
1204C    &   dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
1205C    &   apps=0.
1206C
1207C   limitation des pentes aux poles
1208C     DO ij=1,iip1
1209C        dyq(ij)=appn*dyq(ij)
1210C        dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
1211C     ENDDO
1212C
1213C   test
1214C      DO ij=1,iip1
1215C         dyq(iip1+ij)=0.
1216C         dyq(ip1jm+ij-iip1)=0.
1217C      ENDDO
1218C      DO ij=1,ip1jmp1
1219C         dyq(ij)=dyq(ij)*cos(rlatu((ij-1)/iip1+1))
1220C      ENDDO
1221C
1222C changement 10 07 96
1223C     IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
1224C    &   THEN
1225C        DO ij=1,iip1
1226C           dyqmax(ij)=0.
1227C        ENDDO
1228C     ELSE
1229C        DO ij=1,iip1
1230C           dyqmax(ij)=pente_max*abs(dyqv(ij))
1231C        ENDDO
1232C     ENDIF
1233C
1234C     IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
1235C    & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
1236C    &THEN
1237C        DO ij=ip1jm+1,ip1jmp1
1238C           dyqmax(ij)=0.
1239C        ENDDO
1240C     ELSE
1241C        DO ij=ip1jm+1,ip1jmp1
1242C           dyqmax(ij)=pente_max*abs(dyqv(ij-iip1))
1243C        ENDDO
1244C     ENDIF
1245C   fin changement 10 07 96
1246CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1247
1248c   calcul des pentes limitees
1249      ijb=ij_begin-iip1
1250      ije=ij_end+iip1
1251      if (pole_nord) ijb=ij_begin+iip1
1252      if (pole_sud)  ije=ij_end-iip1
1253
1254      DO ij=ijb,ije
1255         IF(dyqv(ij)*dyqv(ij-iip1).gt.0.) THEN
1256            dyq(ij,l)=sign(min(abs(dyq(ij,l)),dyqmax(ij)),dyq(ij,l))
1257         ELSE
1258            dyq(ij,l)=0.
1259         ENDIF
1260      ENDDO
1261
1262      ENDDO
1263c$OMP END DO NOWAIT
1264
1265      ijb=ij_begin-iip1
1266      ije=ij_end
1267      if (pole_nord) ijb=ij_begin
1268      if (pole_sud)  ije=ij_end-iip1
1269
1270c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1271      DO l=1,llm
1272       DO ij=ijb,ije
1273          IF(masse_adv_v(ij,l).gt.0) THEN
1274              qbyv(ij,l)=q(ij+iip1,l)+dyq(ij+iip1,l)*
1275     ,                   0.5*(1.-masse_adv_v(ij,l)/masse(ij+iip1,l))
1276          ELSE
1277              qbyv(ij,l)=q(ij,l)-dyq(ij,l)*
1278     ,                   0.5*(1.+masse_adv_v(ij,l)/masse(ij,l))
1279          ENDIF
1280          qbyv(ij,l)=masse_adv_v(ij,l)*qbyv(ij,l)
1281       ENDDO
1282      ENDDO
1283c$OMP END DO NOWAIT
1284     
1285      ijb=ij_begin
1286      ije=ij_end
1287      if (pole_nord) ijb=ij_begin+iip1
1288      if (pole_sud)  ije=ij_end-iip1
1289     
1290c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
1291      DO l=1,llm
1292         DO ij=ijb,ije
1293            newmasse=masse(ij,l)
1294     &      +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
1295     
1296            q(ij,l)=(q(ij,l)*masse(ij,l)+qbyv(ij,l)-qbyv(ij-iip1,l))
1297     &         /newmasse
1298            masse(ij,l)=newmasse
1299         ENDDO
1300c.-. ancienne version
1301c        convpn=SSUM(iim,qbyv(1,l),1)/apoln
1302c        convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln
1303         if (pole_nord) then
1304           convpn=SSUM(iim,qbyv(1,l),1)
1305           convmpn=ssum(iim,masse_adv_v(1,l),1)
1306           massepn=ssum(iim,masse(1,l),1)
1307           qpn=0.
1308           do ij=1,iim
1309              qpn=qpn+masse(ij,l)*q(ij,l)
1310           enddo
1311           qpn=(qpn+convpn)/(massepn+convmpn)
1312           do ij=1,iip1
1313              q(ij,l)=qpn
1314           enddo
1315         endif
1316         
1317c        convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)/apols
1318c        convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)/apols
1319         
1320         if (pole_sud) then
1321         
1322           convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
1323           convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
1324           masseps=ssum(iim, masse(ip1jm+1,l),1)
1325           qps=0.
1326           do ij = ip1jm+1,ip1jmp1-1
1327              qps=qps+masse(ij,l)*q(ij,l)
1328           enddo
1329           qps=(qps+convps)/(masseps+convmps)
1330           do ij=ip1jm+1,ip1jmp1
1331              q(ij,l)=qps
1332           enddo
1333         endif
1334c.-. fin ancienne version
1335
1336c._. nouvelle version
1337c        convpn=SSUM(iim,qbyv(1,l),1)
1338c        convmpn=ssum(iim,masse_adv_v(1,l),1)
1339c        oldmasse=ssum(iim,masse(1,l),1)
1340c        newmasse=oldmasse+convmpn
1341c        newq=(q(1,l)*oldmasse+convpn)/newmasse
1342c        newmasse=newmasse/apoln
1343c        DO ij = 1,iip1
1344c           q(ij,l)=newq
1345c           masse(ij,l)=newmasse*aire(ij)
1346c        ENDDO
1347c        convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
1348c        convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
1349c        oldmasse=ssum(iim,masse(ip1jm-iim,l),1)
1350c        newmasse=oldmasse+convmps
1351c        newq=(q(ip1jmp1,l)*oldmasse+convps)/newmasse
1352c        newmasse=newmasse/apols
1353c        DO ij = ip1jm+1,ip1jmp1
1354c           q(ij,l)=newq
1355c           masse(ij,l)=newmasse*aire(ij)
1356c        ENDDO
1357c._. fin nouvelle version
1358      ENDDO
1359c$OMP END DO NOWAIT
1360
1361      RETURN
1362      END
1363     
1364     
1365     
1366      SUBROUTINE vlz_p(q,pente_max,masse,w,ijb_x,ije_x)
1367c
1368c     Auteurs:   P.Le Van, F.Hourdin, F.Forget 
1369c
1370c    ********************************************************************
1371c     Shema  d'advection " pseudo amont " .
1372c    ********************************************************************
1373c    q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....
1374c     dq               sont des arguments de sortie pour le s-pg ....
1375c
1376c
1377c   --------------------------------------------------------------------
1378      USE Parallel_lmdz
1379      IMPLICIT NONE
1380c
1381!-----------------------------------------------------------------------
1382!   INCLUDE 'dimensions.h'
1383!
1384!   dimensions.h contient les dimensions du modele
1385!   ndm est tel que iim=2**ndm
1386!-----------------------------------------------------------------------
1387
1388      INTEGER iim,jjm,llm,ndm
1389
1390      PARAMETER (iim= 128,jjm=96,llm=64,ndm=1)
1391
1392!-----------------------------------------------------------------------
1393!
1394! $Header$
1395!
1396!
1397!  ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
1398!                 veillez  n'utiliser que des ! pour les commentaires
1399!                 et  bien positionner les & des lignes de continuation
1400!                 (les placer en colonne 6 et en colonne 73)
1401!
1402!
1403!-----------------------------------------------------------------------
1404!   INCLUDE 'paramet.h'
1405
1406      INTEGER  iip1,iip2,iip3,jjp1,llmp1,llmp2,llmm1
1407      INTEGER  kftd,ip1jm,ip1jmp1,ip1jmi1,ijp1llm
1408      INTEGER  ijmllm,mvar
1409      INTEGER jcfil,jcfllm
1410
1411      PARAMETER( iip1= iim+1,iip2=iim+2,iip3=iim+3                       &
1412     &    ,jjp1=jjm+1-1/jjm)
1413      PARAMETER( llmp1 = llm+1,  llmp2 = llm+2, llmm1 = llm-1 )
1414      PARAMETER( kftd  = iim/2 -ndm )
1415      PARAMETER( ip1jm  = iip1*jjm,  ip1jmp1= iip1*jjp1 )
1416      PARAMETER( ip1jmi1= ip1jm - iip1 )
1417      PARAMETER( ijp1llm= ip1jmp1 * llm, ijmllm= ip1jm * llm )
1418      PARAMETER( mvar= ip1jmp1*( 2*llm+1) + ijmllm )
1419      PARAMETER( jcfil=jjm/2+5, jcfllm=jcfil*llm )
1420
1421!-----------------------------------------------------------------------
1422!
1423! $Id: logic.h 1520 2011-05-23 11:37:09Z emillour $
1424!
1425!
1426! NB: keep items of different kinds in seperate common blocs to avoid
1427!     "misaligned commons" issues
1428!-----------------------------------------------------------------------
1429! INCLUDE 'logic.h'
1430
1431      COMMON/logicl/ purmats,forward,leapf,apphys,                      &
1432     &  statcl,conser,apdiss,apdelq,saison,ecripar,fxyhypb,ysinus       &
1433     &  ,read_start,ok_guide,ok_strato,tidal,ok_gradsfile               &
1434     &  ,ok_limit,ok_etat0,hybrid                                       &
1435     &  ,moyzon_mu,moyzon_ch
1436
1437      COMMON/logici/ iflag_phys,iflag_trac
1438     
1439      LOGICAL purmats,forward,leapf,apphys,statcl,conser,               &
1440     & apdiss,apdelq,saison,ecripar,fxyhypb,ysinus                      &
1441     &  ,read_start,ok_guide,ok_strato,tidal,ok_gradsfile               &
1442     &  ,ok_limit,ok_etat0
1443      logical hybrid ! vertical coordinate is hybrid if true (sigma otherwise)
1444                     ! (only used if disvert_type==2)
1445      logical moyzon_mu,moyzon_ch ! used for zonal averages in Titan
1446
1447      integer iflag_phys,iflag_trac
1448!$OMP THREADPRIVATE(/logicl/)
1449!$OMP THREADPRIVATE(/logici/)
1450!-----------------------------------------------------------------------
1451!
1452! $Id: comvert.h 1654 2012-09-24 15:07:18Z aslmd $
1453!
1454!-----------------------------------------------------------------------
1455!   INCLUDE 'comvert.h'
1456
1457      COMMON/comvertr/ap(llm+1),bp(llm+1),presnivs(llm),dpres(llm),     &
1458     &               pa,preff,nivsigs(llm),nivsig(llm+1),               &
1459     &               aps(llm),bps(llm),scaleheight,pseudoalt(llm)
1460
1461      common/comverti/disvert_type, pressure_exner
1462
1463      real ap     ! hybrid pressure contribution at interlayers
1464      real bp     ! hybrid sigma contribution at interlayer
1465      real presnivs ! (reference) pressure at mid-layers
1466      real dpres
1467      real pa     ! reference pressure (Pa) at which hybrid coordinates
1468                  ! become purely pressure
1469      real preff  ! reference surface pressure (Pa)
1470      real nivsigs
1471      real nivsig
1472      real aps    ! hybrid pressure contribution at mid-layers
1473      real bps    ! hybrid sigma contribution at mid-layers
1474      real scaleheight ! atmospheric (reference) scale height (km)
1475      real pseudoalt ! pseudo-altitude of model levels (km), based on presnivs(),
1476                     ! preff and scaleheight
1477
1478      integer disvert_type ! type of vertical discretization:
1479                           ! 1: Earth (default for planet_type==earth),
1480                           !     automatic generation
1481                           ! 2: Planets (default for planet_type!=earth),
1482                           !     using 'z2sig.def' (or 'esasig.def) file
1483
1484      logical pressure_exner
1485!     compute pressure inside layers using Exner function, else use mean
1486!     of pressure values at interfaces
1487
1488 !-----------------------------------------------------------------------
1489!
1490! $Id: comconst.h 1437 2010-09-30 08:29:10Z emillour $
1491!
1492!-----------------------------------------------------------------------
1493! INCLUDE comconst.h
1494
1495      COMMON/comconsti/im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl,          &
1496     &                 iflag_top_bound,mode_top_bound
1497      COMMON/comconstr/dtvr,daysec,                                     &
1498     & pi,dtphys,dtdiss,rad,r,kappa,cotot,unsim,g,omeg                  &
1499     & ,dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta        &
1500     & ,dissip_pupstart  ,tau_top_bound,                                &
1501     & daylen,molmass, ihf
1502      COMMON/cpdetvenus/cpp,nu_venus,t0_venus
1503
1504      INTEGER im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl
1505      REAL dtvr ! dynamical time step (in s)
1506      REAL daysec !length (in s) of a standard day
1507      REAL pi    ! something like 3.14159....
1508      REAL dtphys ! (s) time step for the physics
1509      REAL dtdiss ! (s) time step for the dissipation
1510      REAL rad ! (m) radius of the planet
1511      REAL r ! Reduced Gas constant r=R/mu
1512             ! with R=8.31.. J.K-1.mol-1, mu: mol mass of atmosphere (kg/mol)
1513      REAL cpp   ! Cp
1514      REAL kappa ! kappa=R/Cp
1515      REAL cotot
1516      REAL unsim ! = 1./iim
1517      REAL g ! (m/s2) gravity
1518      REAL omeg ! (rad/s) rotation rate of the planet
1519! Dissipation factors, for Earth model:
1520      REAL dissip_factz,dissip_zref !dissip_deltaz
1521! Dissipation factors, for other planets:
1522      REAL dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta
1523      REAL dissip_pupstart
1524      INTEGER iflag_top_bound,mode_top_bound
1525      REAL tau_top_bound
1526      REAL daylen ! length of solar day, in 'standard' day length
1527      REAL molmass ! (g/mol) molar mass of the atmosphere
1528
1529      REAL nu_venus,t0_venus ! coeffs needed for Cp(T), Venus atmosphere
1530      REAL ihf  ! (W/m2) intrinsic heat flux for giant planets
1531
1532
1533!-----------------------------------------------------------------------
1534c
1535c
1536c   Arguments:
1537c   ----------
1538      REAL masse(ip1jmp1,llm),pente_max
1539      REAL q(ip1jmp1,llm)
1540      REAL w(ip1jmp1,llm+1)
1541c
1542c      Local 
1543c   ---------
1544c
1545      INTEGER i,ij,l,j,ii
1546c
1547      REAL,SAVE :: wq(ip1jmp1,llm+1)
1548      REAL newmasse
1549
1550      REAL,SAVE :: dzq(ip1jmp1,llm),dzqw(ip1jmp1,llm),adzqw(ip1jmp1,llm)
1551      REAL dzqmax
1552      REAL sigw
1553
1554      LOGICAL testcpu
1555      SAVE testcpu
1556c$OMP THREADPRIVATE(testcpu)
1557      REAL temps0,temps1,temps2,temps3,temps4,temps5,second
1558      SAVE temps0,temps1,temps2,temps3,temps4,temps5
1559c$OMP THREADPRIVATE(temps0,temps1,temps2,temps3,temps4,temps5)
1560
1561      REAL      SSUM
1562      EXTERNAL  SSUM
1563
1564      DATA testcpu/.false./
1565      DATA temps0,temps1,temps2,temps3,temps4,temps5/0.,0.,0.,0.,0.,0./
1566      INTEGER ijb,ije,ijb_x,ije_x
1567c    On oriente tout dans le sens de la pression c'est a dire dans le
1568c    sens de W
1569
1570
1571      ijb=ijb_x
1572      ije=ije_x
1573
1574c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
1575      DO l=2,llm
1576         DO ij=ijb,ije
1577            dzqw(ij,l)=q(ij,l-1)-q(ij,l)
1578            adzqw(ij,l)=abs(dzqw(ij,l))
1579         ENDDO
1580      ENDDO
1581c$OMP END DO
1582
1583c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1584      DO l=2,llm-1
1585         DO ij=ijb,ije
1586            IF(dzqw(ij,l)*dzqw(ij,l+1).gt.0.) THEN
1587                dzq(ij,l)=0.5*(dzqw(ij,l)+dzqw(ij,l+1))
1588            ELSE
1589                dzq(ij,l)=0.
1590            ENDIF
1591            dzqmax=pente_max*min(adzqw(ij,l),adzqw(ij,l+1))
1592            dzq(ij,l)=sign(min(abs(dzq(ij,l)),dzqmax),dzq(ij,l))
1593         ENDDO
1594      ENDDO
1595c$OMP END DO NOWAIT
1596
1597c$OMP MASTER
1598      DO ij=ijb,ije
1599         dzq(ij,1)=0.
1600         dzq(ij,llm)=0.
1601      ENDDO
1602c$OMP END MASTER
1603c$OMP BARRIER
1604c ---------------------------------------------------------------
1605c   .... calcul des termes d'advection verticale  .......
1606c ---------------------------------------------------------------
1607
1608c calcul de  - d( q   * w )/ d(sigma)    qu'on ajoute a  dq pour calculer dq
1609
1610c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1611       DO l = 1,llm-1
1612         do  ij = ijb,ije
1613          IF(w(ij,l+1).gt.0.) THEN
1614             sigw=w(ij,l+1)/masse(ij,l+1)
1615             wq(ij,l+1)=w(ij,l+1)*(q(ij,l+1)+0.5*(1.-sigw)*dzq(ij,l+1))
1616          ELSE
1617             sigw=w(ij,l+1)/masse(ij,l)
1618             wq(ij,l+1)=w(ij,l+1)*(q(ij,l)-0.5*(1.+sigw)*dzq(ij,l))
1619          ENDIF
1620         ENDDO
1621       ENDDO
1622c$OMP END DO NOWAIT
1623
1624c$OMP MASTER
1625       DO ij=ijb,ije
1626          wq(ij,llm+1)=0.
1627          wq(ij,1)=0.
1628       ENDDO
1629c$OMP END MASTER
1630c$OMP BARRIER
1631
1632c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1633      DO l=1,llm
1634         DO ij=ijb,ije
1635            newmasse=masse(ij,l)+w(ij,l+1)-w(ij,l)
1636            q(ij,l)=(q(ij,l)*masse(ij,l)+wq(ij,l+1)-wq(ij,l))
1637     &         /newmasse
1638            masse(ij,l)=newmasse
1639         ENDDO
1640      ENDDO
1641c$OMP END DO NOWAIT
1642
1643
1644      RETURN
1645      END
1646c      SUBROUTINE minmaxq(zq,qmin,qmax,comment)
1647c
1648c#include "dimensions.h"
1649c#include "paramet.h"
1650
1651c      CHARACTER*(*) comment
1652c      real qmin,qmax
1653c      real zq(ip1jmp1,llm)
1654
1655c      INTEGER jadrs(ip1jmp1), jbad, k, i
1656
1657
1658c      DO k = 1, llm
1659c         jbad = 0
1660c         DO i = 1, ip1jmp1
1661c         IF (zq(i,k).GT.qmax .OR. zq(i,k).LT.qmin) THEN
1662c            jbad = jbad + 1
1663c            jadrs(jbad) = i
1664c         ENDIF
1665c         ENDDO
1666c         IF (jbad.GT.0) THEN
1667c         PRINT*, comment
1668c         DO i = 1, jbad
1669cc            PRINT*, "i,k,zq=", jadrs(i),k,zq(jadrs(i),k)
1670c         ENDDO
1671c         ENDIF
1672c      ENDDO
1673
1674c      return
1675c      end
1676
1677
1678      subroutine minmaxq_p(zq,qmin,qmax,comment)
1679
1680!-----------------------------------------------------------------------
1681!   INCLUDE 'dimensions.h'
1682!
1683!   dimensions.h contient les dimensions du modele
1684!   ndm est tel que iim=2**ndm
1685!-----------------------------------------------------------------------
1686
1687      INTEGER iim,jjm,llm,ndm
1688
1689      PARAMETER (iim= 128,jjm=96,llm=64,ndm=1)
1690
1691!-----------------------------------------------------------------------
1692!
1693! $Header$
1694!
1695!
1696!  ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
1697!                 veillez  n'utiliser que des ! pour les commentaires
1698!                 et  bien positionner les & des lignes de continuation
1699!                 (les placer en colonne 6 et en colonne 73)
1700!
1701!
1702!-----------------------------------------------------------------------
1703!   INCLUDE 'paramet.h'
1704
1705      INTEGER  iip1,iip2,iip3,jjp1,llmp1,llmp2,llmm1
1706      INTEGER  kftd,ip1jm,ip1jmp1,ip1jmi1,ijp1llm
1707      INTEGER  ijmllm,mvar
1708      INTEGER jcfil,jcfllm
1709
1710      PARAMETER( iip1= iim+1,iip2=iim+2,iip3=iim+3                       &
1711     &    ,jjp1=jjm+1-1/jjm)
1712      PARAMETER( llmp1 = llm+1,  llmp2 = llm+2, llmm1 = llm-1 )
1713      PARAMETER( kftd  = iim/2 -ndm )
1714      PARAMETER( ip1jm  = iip1*jjm,  ip1jmp1= iip1*jjp1 )
1715      PARAMETER( ip1jmi1= ip1jm - iip1 )
1716      PARAMETER( ijp1llm= ip1jmp1 * llm, ijmllm= ip1jm * llm )
1717      PARAMETER( mvar= ip1jmp1*( 2*llm+1) + ijmllm )
1718      PARAMETER( jcfil=jjm/2+5, jcfllm=jcfil*llm )
1719
1720!-----------------------------------------------------------------------
1721
1722      character*20 comment
1723      real qmin,qmax
1724      real zq(ip1jmp1,llm)
1725      real zzq(iip1,jjp1,llm)
1726
1727      integer imin,jmin,lmin,ijlmin
1728      integer imax,jmax,lmax,ijlmax
1729
1730      integer ismin,ismax
1731
1732      return
17339999  format(a20,'  q(',i3,',',i2,',',i2,')=',e12.5,e12.5)
1734      end
1735
1736
1737
Note: See TracBrowser for help on using the repository browser.