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

Last change on this file since 224 was 224, checked in by ymipsl, 10 years ago
File size: 44.5 KB
Line 
1!
2! $Header$
3!
4       SUBROUTINE vlspltqs_p ( q,pente_max,masse,w,pbaru,pbarv,pdt,
5     ,                                  p,pk,teta                 )
6c
7c     Auteurs:   P.Le Van, F.Hourdin, F.Forget, F.Codron 
8c
9c    ********************************************************************
10c          Shema  d'advection " pseudo amont " .
11c      + test sur humidite specifique: Q advecte< Qsat aval
12c                   (F. Codron, 10/99)
13c    ********************************************************************
14c     q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....
15c
16c     pente_max facteur de limitation des pentes: 2 en general
17c                                                0 pour un schema amont
18c     pbaru,pbarv,w flux de masse en u ,v ,w
19c     pdt pas de temps
20c
21c     teta temperature potentielle, p pression aux interfaces,
22c     pk exner au milieu des couches necessaire pour calculer Qsat
23c   --------------------------------------------------------------------
24      USE parallel_lmdz
25      USE mod_hallo
26      USE VAMPIR
27      use cpdet_mod, only: tpot2t_glo_p
28      IMPLICIT NONE
29
30c
31!-----------------------------------------------------------------------
32!   INCLUDE 'dimensions.h'
33!
34!   dimensions.h contient les dimensions du modele
35!   ndm est tel que iim=2**ndm
36!-----------------------------------------------------------------------
37
38      INTEGER iim,jjm,llm,ndm
39
40      PARAMETER (iim= 128,jjm=96,llm=64,ndm=1)
41
42!-----------------------------------------------------------------------
43!
44! $Header$
45!
46!
47!  ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
48!                 veillez  n'utiliser que des ! pour les commentaires
49!                 et  bien positionner les & des lignes de continuation
50!                 (les placer en colonne 6 et en colonne 73)
51!
52!
53!-----------------------------------------------------------------------
54!   INCLUDE 'paramet.h'
55
56      INTEGER  iip1,iip2,iip3,jjp1,llmp1,llmp2,llmm1
57      INTEGER  kftd,ip1jm,ip1jmp1,ip1jmi1,ijp1llm
58      INTEGER  ijmllm,mvar
59      INTEGER jcfil,jcfllm
60
61      PARAMETER( iip1= iim+1,iip2=iim+2,iip3=iim+3                       &
62     &    ,jjp1=jjm+1-1/jjm)
63      PARAMETER( llmp1 = llm+1,  llmp2 = llm+2, llmm1 = llm-1 )
64      PARAMETER( kftd  = iim/2 -ndm )
65      PARAMETER( ip1jm  = iip1*jjm,  ip1jmp1= iip1*jjp1 )
66      PARAMETER( ip1jmi1= ip1jm - iip1 )
67      PARAMETER( ijp1llm= ip1jmp1 * llm, ijmllm= ip1jm * llm )
68      PARAMETER( mvar= ip1jmp1*( 2*llm+1) + ijmllm )
69      PARAMETER( jcfil=jjm/2+5, jcfllm=jcfil*llm )
70
71!-----------------------------------------------------------------------
72!
73! $Id: logic.h 1520 2011-05-23 11:37:09Z emillour $
74!
75!
76! NB: keep items of different kinds in seperate common blocs to avoid
77!     "misaligned commons" issues
78!-----------------------------------------------------------------------
79! INCLUDE 'logic.h'
80
81      COMMON/logicl/ purmats,forward,leapf,apphys,                      &
82     &  statcl,conser,apdiss,apdelq,saison,ecripar,fxyhypb,ysinus       &
83     &  ,read_start,ok_guide,ok_strato,tidal,ok_gradsfile               &
84     &  ,ok_limit,ok_etat0,hybrid                                       &
85     &  ,moyzon_mu,moyzon_ch
86
87      COMMON/logici/ iflag_phys,iflag_trac
88     
89      LOGICAL purmats,forward,leapf,apphys,statcl,conser,               &
90     & apdiss,apdelq,saison,ecripar,fxyhypb,ysinus                      &
91     &  ,read_start,ok_guide,ok_strato,tidal,ok_gradsfile               &
92     &  ,ok_limit,ok_etat0
93      logical hybrid ! vertical coordinate is hybrid if true (sigma otherwise)
94                     ! (only used if disvert_type==2)
95      logical moyzon_mu,moyzon_ch ! used for zonal averages in Titan
96
97      integer iflag_phys,iflag_trac
98!$OMP THREADPRIVATE(/logicl/)
99!$OMP THREADPRIVATE(/logici/)
100!-----------------------------------------------------------------------
101!
102! $Id: comvert.h 1654 2012-09-24 15:07:18Z aslmd $
103!
104!-----------------------------------------------------------------------
105!   INCLUDE 'comvert.h'
106
107      COMMON/comvertr/ap(llm+1),bp(llm+1),presnivs(llm),dpres(llm),     &
108     &               pa,preff,nivsigs(llm),nivsig(llm+1),               &
109     &               aps(llm),bps(llm),scaleheight,pseudoalt(llm)
110
111      common/comverti/disvert_type, pressure_exner
112
113      real ap     ! hybrid pressure contribution at interlayers
114      real bp     ! hybrid sigma contribution at interlayer
115      real presnivs ! (reference) pressure at mid-layers
116      real dpres
117      real pa     ! reference pressure (Pa) at which hybrid coordinates
118                  ! become purely pressure
119      real preff  ! reference surface pressure (Pa)
120      real nivsigs
121      real nivsig
122      real aps    ! hybrid pressure contribution at mid-layers
123      real bps    ! hybrid sigma contribution at mid-layers
124      real scaleheight ! atmospheric (reference) scale height (km)
125      real pseudoalt ! pseudo-altitude of model levels (km), based on presnivs(),
126                     ! preff and scaleheight
127
128      integer disvert_type ! type of vertical discretization:
129                           ! 1: Earth (default for planet_type==earth),
130                           !     automatic generation
131                           ! 2: Planets (default for planet_type!=earth),
132                           !     using 'z2sig.def' (or 'esasig.def) file
133
134      logical pressure_exner
135!     compute pressure inside layers using Exner function, else use mean
136!     of pressure values at interfaces
137
138 !-----------------------------------------------------------------------
139!
140! $Id: comconst.h 1437 2010-09-30 08:29:10Z emillour $
141!
142!-----------------------------------------------------------------------
143! INCLUDE comconst.h
144
145      COMMON/comconsti/im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl,          &
146     &                 iflag_top_bound,mode_top_bound
147      COMMON/comconstr/dtvr,daysec,                                     &
148     & pi,dtphys,dtdiss,rad,r,kappa,cotot,unsim,g,omeg                  &
149     & ,dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta        &
150     & ,dissip_pupstart  ,tau_top_bound,                                &
151     & daylen,molmass, ihf
152      COMMON/cpdetvenus/cpp,nu_venus,t0_venus
153
154      INTEGER im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl
155      REAL dtvr ! dynamical time step (in s)
156      REAL daysec !length (in s) of a standard day
157      REAL pi    ! something like 3.14159....
158      REAL dtphys ! (s) time step for the physics
159      REAL dtdiss ! (s) time step for the dissipation
160      REAL rad ! (m) radius of the planet
161      REAL r ! Reduced Gas constant r=R/mu
162             ! with R=8.31.. J.K-1.mol-1, mu: mol mass of atmosphere (kg/mol)
163      REAL cpp   ! Cp
164      REAL kappa ! kappa=R/Cp
165      REAL cotot
166      REAL unsim ! = 1./iim
167      REAL g ! (m/s2) gravity
168      REAL omeg ! (rad/s) rotation rate of the planet
169! Dissipation factors, for Earth model:
170      REAL dissip_factz,dissip_zref !dissip_deltaz
171! Dissipation factors, for other planets:
172      REAL dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta
173      REAL dissip_pupstart
174      INTEGER iflag_top_bound,mode_top_bound
175      REAL tau_top_bound
176      REAL daylen ! length of solar day, in 'standard' day length
177      REAL molmass ! (g/mol) molar mass of the atmosphere
178
179      REAL nu_venus,t0_venus ! coeffs needed for Cp(T), Venus atmosphere
180      REAL ihf  ! (W/m2) intrinsic heat flux for giant planets
181
182
183!-----------------------------------------------------------------------
184
185c
186c   Arguments:
187c   ----------
188      REAL masse(ip1jmp1,llm),pente_max
189      REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
190      REAL q(ip1jmp1,llm)
191      REAL w(ip1jmp1,llm),pdt
192      REAL p(ip1jmp1,llmp1),teta(ip1jmp1,llm),pk(ip1jmp1,llm)
193c
194c      Local 
195c   ---------
196c
197      INTEGER i,ij,l,j,ii
198c
199      REAL qsat(ip1jmp1,llm)
200      REAL zm(ip1jmp1,llm)
201      REAL mu(ip1jmp1,llm)
202      REAL mv(ip1jm,llm)
203      REAL mw(ip1jmp1,llm+1)
204      REAL zq(ip1jmp1,llm)
205      REAL temps1,temps2,temps3
206      REAL zzpbar, zzw
207      LOGICAL testcpu
208      SAVE testcpu
209      SAVE temps1,temps2,temps3
210
211      REAL qmin,qmax
212      DATA qmin,qmax/0.,1.e33/
213      DATA testcpu/.false./
214      DATA temps1,temps2,temps3/0.,0.,0./
215
216c--pour rapport de melange saturant--
217
218      REAL rtt,retv,r2es,r3les,r3ies,r4les,r4ies,play
219      REAL ptarg,pdelarg,foeew,zdelta
220!      REAL tempe(ip1jmp1)
221! ADAPTATION GCM POUR CP(T)
222      REAL tempe(ip1jmp1,llm)
223
224      INTEGER ijb,ije
225      type(request) :: MyRequest1
226      type(request) :: MyRequest2
227
228c    fonction psat(T)
229
230       FOEEW ( PTARG,PDELARG ) = EXP (
231     *          (R3LES*(1.-PDELARG)+R3IES*PDELARG) * (PTARG-RTT)
232     * / (PTARG-(R4LES*(1.-PDELARG)+R4IES*PDELARG)) )
233
234        r2es  = 380.11733 
235        r3les = 17.269
236        r3ies = 21.875
237        r4les = 35.86
238        r4ies = 7.66
239        retv = 0.6077667
240        rtt  = 273.16
241
242c-- Calcul de Qsat en chaque point
243c-- approximation: au milieu des couches play(l)=(p(l)+p(l+1))/2
244c   pour eviter une exponentielle.
245
246! ADAPTATION GCM POUR CP(T)
247! probablement a revoir...
248!      call tpot2t_p(ip1jmp1,llm,teta,tempe,pk)
249      call tpot2t_glo_p(teta,tempe,pk)
250
251      call SetTag(MyRequest1,100)
252      call SetTag(MyRequest2,101)
253       
254        ijb=ij_begin-iip1
255        ije=ij_end+iip1
256        if (pole_nord) ijb=ij_begin
257        if (pole_sud) ije=ij_end
258       
259       
260        DO l = 1, llm
261!         DO ij = ijb, ije
262!          tempe(ij) = teta(ij,l) * pk(ij,l) /cpp
263!         ENDDO
264         DO ij = ijb, ije
265          zdelta = MAX( 0., SIGN(1., rtt - tempe(ij,l)) )
266          play   = 0.5*(p(ij,l)+p(ij,l+1))
267          qsat(ij,l) = MIN(0.5, r2es* FOEEW(tempe(ij,l),zdelta) / play )
268          qsat(ij,l) = qsat(ij,l) / ( 1. - retv * qsat(ij,l) )
269         ENDDO
270        ENDDO
271
272c      PRINT*,'Debut vlsplt version debug sans vlyqs'
273
274        zzpbar = 0.5 * pdt
275        zzw    = pdt
276
277      ijb=ij_begin
278      ije=ij_end
279      if (pole_nord) ijb=ijb+iip1
280      if (pole_sud)  ije=ije-iip1
281
282     
283      DO l=1,llm
284        DO ij = ijb,ije
285            mu(ij,l)=pbaru(ij,l) * zzpbar
286         ENDDO
287      ENDDO
288     
289      ijb=ij_begin-iip1
290      ije=ij_end
291      if (pole_nord) ijb=ij_begin
292      if (pole_sud)  ije=ij_end-iip1
293
294      DO l=1,llm
295         DO ij=ijb,ije
296            mv(ij,l)=pbarv(ij,l) * zzpbar
297         ENDDO
298      ENDDO
299
300      ijb=ij_begin
301      ije=ij_end
302     
303      DO l=1,llm
304         DO ij=ijb,ije
305            mw(ij,l)=w(ij,l) * zzw
306         ENDDO
307      ENDDO
308
309      DO ij=ijb,ije
310         mw(ij,llm+1)=0.
311      ENDDO
312
313c      CALL SCOPY(ijp1llm,q,1,zq,1)
314c      CALL SCOPY(ijp1llm,masse,1,zm,1)
315
316       ijb=ij_begin
317       ije=ij_end
318       zq(ijb:ije,1:llm)=q(ijb:ije,1:llm)
319       zm(ijb:ije,1:llm)=masse(ijb:ije,1:llm)
320
321
322      call vlxqs_p(zq,pente_max,zm,mu,qsat,ij_begin,ij_begin+2*iip1-1)
323      call vlxqs_p(zq,pente_max,zm,mu,qsat,ij_end-2*iip1+1,ij_end)
324 
325      call VTb(VTHallo)
326      call Register_Hallo(zq,ip1jmp1,llm,2,2,2,2,MyRequest1)
327      call Register_Hallo(zm,ip1jmp1,llm,1,1,1,1,MyRequest1)
328      call SendRequest(MyRequest1)
329      call VTe(VTHallo)
330
331      call vlxqs_p(zq,pente_max,zm,mu,qsat,
332     .             ij_begin+2*iip1,ij_end-2*iip1)
333
334      call VTb(VTHallo)
335      call WaitRecvRequest(MyRequest1)
336      call VTe(VTHallo)
337
338      call vlyqs_p(zq,pente_max,zm,mv,qsat)
339
340      call vlz_p(zq,pente_max,zm,mw,ij_begin,ij_begin+2*iip1-1)
341      call vlz_p(zq,pente_max,zm,mw,ij_end-2*iip1+1,ij_end)
342
343      call VTb(VTHallo)
344      call Register_Hallo(zq,ip1jmp1,llm,2,2,2,2,MyRequest2)
345      call Register_Hallo(zm,ip1jmp1,llm,1,1,1,1,MyRequest2)
346      call SendRequest(MyRequest2)
347      call VTe(VTHallo)
348
349      call vlz_p(zq,pente_max,zm,mw,ij_begin+2*iip1,ij_end-2*iip1)
350
351      call VTb(VTHallo)
352      call WaitRecvRequest(MyRequest2)
353      call VTe(VTHallo)
354     
355      call vlyqs_p(zq,pente_max,zm,mv,qsat)
356
357
358      call vlxqs_p(zq,pente_max,zm,mu,qsat,ij_begin,ij_end)
359
360
361      ijb=ij_begin
362      ije=ij_end
363
364      DO l=1,llm
365         DO ij=ijb,ije
366           q(ij,l)=zq(ij,l)
367         ENDDO
368      ENDDO
369     
370      DO l=1,llm
371         DO ij=ijb,ije-iip1+1,iip1
372            q(ij+iim,l)=q(ij,l)
373         ENDDO
374      ENDDO
375
376      call WaitSendRequest(MyRequest1) 
377      call WaitSendRequest(MyRequest2)
378
379      RETURN
380      END
381      SUBROUTINE vlxqs_p(q,pente_max,masse,u_m,qsat,ijb_x,ije_x)
382c
383c     Auteurs:   P.Le Van, F.Hourdin, F.Forget 
384c
385c    ********************************************************************
386c     Shema  d'advection " pseudo amont " .
387c    ********************************************************************
388c
389c   --------------------------------------------------------------------
390      USE parallel_lmdz
391      IMPLICIT NONE
392c
393!-----------------------------------------------------------------------
394!   INCLUDE 'dimensions.h'
395!
396!   dimensions.h contient les dimensions du modele
397!   ndm est tel que iim=2**ndm
398!-----------------------------------------------------------------------
399
400      INTEGER iim,jjm,llm,ndm
401
402      PARAMETER (iim= 128,jjm=96,llm=64,ndm=1)
403
404!-----------------------------------------------------------------------
405!
406! $Header$
407!
408!
409!  ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
410!                 veillez  n'utiliser que des ! pour les commentaires
411!                 et  bien positionner les & des lignes de continuation
412!                 (les placer en colonne 6 et en colonne 73)
413!
414!
415!-----------------------------------------------------------------------
416!   INCLUDE 'paramet.h'
417
418      INTEGER  iip1,iip2,iip3,jjp1,llmp1,llmp2,llmm1
419      INTEGER  kftd,ip1jm,ip1jmp1,ip1jmi1,ijp1llm
420      INTEGER  ijmllm,mvar
421      INTEGER jcfil,jcfllm
422
423      PARAMETER( iip1= iim+1,iip2=iim+2,iip3=iim+3                       &
424     &    ,jjp1=jjm+1-1/jjm)
425      PARAMETER( llmp1 = llm+1,  llmp2 = llm+2, llmm1 = llm-1 )
426      PARAMETER( kftd  = iim/2 -ndm )
427      PARAMETER( ip1jm  = iip1*jjm,  ip1jmp1= iip1*jjp1 )
428      PARAMETER( ip1jmi1= ip1jm - iip1 )
429      PARAMETER( ijp1llm= ip1jmp1 * llm, ijmllm= ip1jm * llm )
430      PARAMETER( mvar= ip1jmp1*( 2*llm+1) + ijmllm )
431      PARAMETER( jcfil=jjm/2+5, jcfllm=jcfil*llm )
432
433!-----------------------------------------------------------------------
434!
435! $Id: logic.h 1520 2011-05-23 11:37:09Z emillour $
436!
437!
438! NB: keep items of different kinds in seperate common blocs to avoid
439!     "misaligned commons" issues
440!-----------------------------------------------------------------------
441! INCLUDE 'logic.h'
442
443      COMMON/logicl/ purmats,forward,leapf,apphys,                      &
444     &  statcl,conser,apdiss,apdelq,saison,ecripar,fxyhypb,ysinus       &
445     &  ,read_start,ok_guide,ok_strato,tidal,ok_gradsfile               &
446     &  ,ok_limit,ok_etat0,hybrid                                       &
447     &  ,moyzon_mu,moyzon_ch
448
449      COMMON/logici/ iflag_phys,iflag_trac
450     
451      LOGICAL purmats,forward,leapf,apphys,statcl,conser,               &
452     & apdiss,apdelq,saison,ecripar,fxyhypb,ysinus                      &
453     &  ,read_start,ok_guide,ok_strato,tidal,ok_gradsfile               &
454     &  ,ok_limit,ok_etat0
455      logical hybrid ! vertical coordinate is hybrid if true (sigma otherwise)
456                     ! (only used if disvert_type==2)
457      logical moyzon_mu,moyzon_ch ! used for zonal averages in Titan
458
459      integer iflag_phys,iflag_trac
460!$OMP THREADPRIVATE(/logicl/)
461!$OMP THREADPRIVATE(/logici/)
462!-----------------------------------------------------------------------
463!
464! $Id: comvert.h 1654 2012-09-24 15:07:18Z aslmd $
465!
466!-----------------------------------------------------------------------
467!   INCLUDE 'comvert.h'
468
469      COMMON/comvertr/ap(llm+1),bp(llm+1),presnivs(llm),dpres(llm),     &
470     &               pa,preff,nivsigs(llm),nivsig(llm+1),               &
471     &               aps(llm),bps(llm),scaleheight,pseudoalt(llm)
472
473      common/comverti/disvert_type, pressure_exner
474
475      real ap     ! hybrid pressure contribution at interlayers
476      real bp     ! hybrid sigma contribution at interlayer
477      real presnivs ! (reference) pressure at mid-layers
478      real dpres
479      real pa     ! reference pressure (Pa) at which hybrid coordinates
480                  ! become purely pressure
481      real preff  ! reference surface pressure (Pa)
482      real nivsigs
483      real nivsig
484      real aps    ! hybrid pressure contribution at mid-layers
485      real bps    ! hybrid sigma contribution at mid-layers
486      real scaleheight ! atmospheric (reference) scale height (km)
487      real pseudoalt ! pseudo-altitude of model levels (km), based on presnivs(),
488                     ! preff and scaleheight
489
490      integer disvert_type ! type of vertical discretization:
491                           ! 1: Earth (default for planet_type==earth),
492                           !     automatic generation
493                           ! 2: Planets (default for planet_type!=earth),
494                           !     using 'z2sig.def' (or 'esasig.def) file
495
496      logical pressure_exner
497!     compute pressure inside layers using Exner function, else use mean
498!     of pressure values at interfaces
499
500 !-----------------------------------------------------------------------
501!
502! $Id: comconst.h 1437 2010-09-30 08:29:10Z emillour $
503!
504!-----------------------------------------------------------------------
505! INCLUDE comconst.h
506
507      COMMON/comconsti/im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl,          &
508     &                 iflag_top_bound,mode_top_bound
509      COMMON/comconstr/dtvr,daysec,                                     &
510     & pi,dtphys,dtdiss,rad,r,kappa,cotot,unsim,g,omeg                  &
511     & ,dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta        &
512     & ,dissip_pupstart  ,tau_top_bound,                                &
513     & daylen,molmass, ihf
514      COMMON/cpdetvenus/cpp,nu_venus,t0_venus
515
516      INTEGER im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl
517      REAL dtvr ! dynamical time step (in s)
518      REAL daysec !length (in s) of a standard day
519      REAL pi    ! something like 3.14159....
520      REAL dtphys ! (s) time step for the physics
521      REAL dtdiss ! (s) time step for the dissipation
522      REAL rad ! (m) radius of the planet
523      REAL r ! Reduced Gas constant r=R/mu
524             ! with R=8.31.. J.K-1.mol-1, mu: mol mass of atmosphere (kg/mol)
525      REAL cpp   ! Cp
526      REAL kappa ! kappa=R/Cp
527      REAL cotot
528      REAL unsim ! = 1./iim
529      REAL g ! (m/s2) gravity
530      REAL omeg ! (rad/s) rotation rate of the planet
531! Dissipation factors, for Earth model:
532      REAL dissip_factz,dissip_zref !dissip_deltaz
533! Dissipation factors, for other planets:
534      REAL dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta
535      REAL dissip_pupstart
536      INTEGER iflag_top_bound,mode_top_bound
537      REAL tau_top_bound
538      REAL daylen ! length of solar day, in 'standard' day length
539      REAL molmass ! (g/mol) molar mass of the atmosphere
540
541      REAL nu_venus,t0_venus ! coeffs needed for Cp(T), Venus atmosphere
542      REAL ihf  ! (W/m2) intrinsic heat flux for giant planets
543
544
545!-----------------------------------------------------------------------
546c
547c
548c   Arguments:
549c   ----------
550      REAL masse(ip1jmp1,llm),pente_max
551      REAL u_m( ip1jmp1,llm )
552      REAL q(ip1jmp1,llm)
553      REAL qsat(ip1jmp1,llm)
554c
555c      Local 
556c   ---------
557c
558      INTEGER ij,l,j,i,iju,ijq,indu(ip1jmp1),niju
559      INTEGER n0,iadvplus(ip1jmp1,llm),nl(llm)
560c
561      REAL new_m,zu_m,zdum(ip1jmp1,llm)
562      REAL dxq(ip1jmp1,llm),dxqu(ip1jmp1)
563      REAL zz(ip1jmp1)
564      REAL adxqu(ip1jmp1),dxqmax(ip1jmp1,llm)
565      REAL u_mq(ip1jmp1,llm)
566
567      REAL      SSUM
568
569
570      INTEGER ijb,ije,ijb_x,ije_x
571     
572
573c   calcul de la pente a droite et a gauche de la maille
574
575c      ijb=ij_begin
576c      ije=ij_end
577
578      ijb=ijb_x
579      ije=ije_x
580       
581      if (pole_nord.and.ijb==1) ijb=ijb+iip1
582      if (pole_sud.and.ije==ip1jmp1)  ije=ije-iip1
583     
584      IF (pente_max.gt.-1.e-5) THEN
585c     IF (pente_max.gt.10) THEN
586
587c   calcul des pentes avec limitation, Van Leer scheme I:
588c   -----------------------------------------------------
589
590c   calcul de la pente aux points u
591
592c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 
593         DO l = 1, llm
594            DO ij=ijb,ije-1
595               dxqu(ij)=q(ij+1,l)-q(ij,l)
596c              IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
597c              sigu(ij)=u_m(ij,l)/masse(ij,l)
598            ENDDO
599            DO ij=ijb+iip1-1,ije,iip1
600               dxqu(ij)=dxqu(ij-iim)
601c              sigu(ij)=sigu(ij-iim)
602            ENDDO
603
604            DO ij=ijb,ije
605               adxqu(ij)=abs(dxqu(ij))
606            ENDDO
607
608c   calcul de la pente maximum dans la maille en valeur absolue
609
610            DO ij=ijb+1,ije
611               dxqmax(ij,l)=pente_max*
612     ,      min(adxqu(ij-1),adxqu(ij))
613c limitation subtile
614c    ,      min(adxqu(ij-1)/sigu(ij-1),adxqu(ij)/(1.-sigu(ij)))
615         
616
617            ENDDO
618
619            DO ij=ijb+iip1-1,ije,iip1
620               dxqmax(ij-iim,l)=dxqmax(ij,l)
621            ENDDO
622
623            DO ij=ijb+1,ije
624               IF(dxqu(ij-1)*dxqu(ij).gt.0) THEN
625                  dxq(ij,l)=dxqu(ij-1)+dxqu(ij)
626               ELSE
627c   extremum local
628                  dxq(ij,l)=0.
629               ENDIF
630               dxq(ij,l)=0.5*dxq(ij,l)
631               dxq(ij,l)=
632     ,         sign(min(abs(dxq(ij,l)),dxqmax(ij,l)),dxq(ij,l))
633            ENDDO
634
635         ENDDO ! l=1,llm
636c$OMP END DO NOWAIT
637
638      ELSE ! (pente_max.lt.-1.e-5)
639
640c   Pentes produits:
641c   ----------------
642c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 
643         DO l = 1, llm
644            DO ij=ijb,ije-1
645               dxqu(ij)=q(ij+1,l)-q(ij,l)
646            ENDDO
647            DO ij=ijb+iip1-1,ije,iip1
648               dxqu(ij)=dxqu(ij-iim)
649            ENDDO
650
651            DO ij=ijb+1,ije
652               zz(ij)=dxqu(ij-1)*dxqu(ij)
653               zz(ij)=zz(ij)+zz(ij)
654               IF(zz(ij).gt.0) THEN
655                  dxq(ij,l)=zz(ij)/(dxqu(ij-1)+dxqu(ij))
656               ELSE
657c   extremum local
658                  dxq(ij,l)=0.
659               ENDIF
660            ENDDO
661
662         ENDDO
663c$OMP END DO NOWAIT
664      ENDIF ! (pente_max.lt.-1.e-5)
665
666c   bouclage de la pente en iip1:
667c   -----------------------------
668c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
669      DO l=1,llm
670         DO ij=ijb+iip1-1,ije,iip1
671            dxq(ij-iim,l)=dxq(ij,l)
672         ENDDO
673
674         DO ij=ijb,ije
675            iadvplus(ij,l)=0
676         ENDDO
677
678      ENDDO
679c$OMP END DO NOWAIT
680     
681      if (pole_nord) THEN
682c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
683        DO l=1,llm     
684          iadvplus(1:iip1,l)=0
685        ENDDO
686c$OMP END DO NOWAIT
687      endif
688     
689      if (pole_sud)  THEN
690c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
691        DO l=1,llm 
692          iadvplus(ip1jm+1:ip1jmp1,l)=0
693        ENDDO
694c$OMP END DO NOWAIT
695      endif
696       
697c   calcul des flux a gauche et a droite
698
699c   on cumule le flux correspondant a toutes les mailles dont la masse
700c   au travers de la paroi pENDant le pas de temps.
701c   le rapport de melange de l'air advecte est min(q_vanleer, Qsat_downwind)
702c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
703      DO l=1,llm
704       DO ij=ijb,ije-1
705          IF (u_m(ij,l).gt.0.) THEN
706             zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l)
707             u_mq(ij,l)=u_m(ij,l)*
708     $         min(q(ij,l)+0.5*zdum(ij,l)*dxq(ij,l),qsat(ij+1,l))
709          ELSE
710             zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l)
711             u_mq(ij,l)=u_m(ij,l)*
712     $         min(q(ij+1,l)-0.5*zdum(ij,l)*dxq(ij+1,l),qsat(ij,l))
713          ENDIF
714       ENDDO
715      ENDDO
716c$OMP END DO NOWAIT
717
718
719c   detection des points ou on advecte plus que la masse de la
720c   maille
721c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
722      DO l=1,llm
723         DO ij=ijb,ije-1
724            IF(zdum(ij,l).lt.0) THEN
725               iadvplus(ij,l)=1
726               u_mq(ij,l)=0.
727            ENDIF
728         ENDDO
729      ENDDO
730c$OMP END DO NOWAIT
731
732c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
733      DO l=1,llm
734       DO ij=ijb+iip1-1,ije,iip1
735          iadvplus(ij,l)=iadvplus(ij-iim,l)
736       ENDDO
737      ENDDO
738c$OMP END DO NOWAIT
739
740
741
742c   traitement special pour le cas ou on advecte en longitude plus que le
743c   contenu de la maille.
744c   cette partie est mal vectorisee.
745
746c   pas d'influence de la pression saturante (pour l'instant)
747
748c  calcul du nombre de maille sur lequel on advecte plus que la maille.
749
750      n0=0
751c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
752      DO l=1,llm
753         nl(l)=0
754         DO ij=ijb,ije
755            nl(l)=nl(l)+iadvplus(ij,l)
756         ENDDO
757         n0=n0+nl(l)
758      ENDDO
759c$OMP END DO NOWAIT
760
761cym ATTENTION ICI en OpenMP reduction pas forcement nécessaire
762cym      IF(n0.gt.1) THEN
763cym        IF(n0.gt.0) THEN
764ccc      PRINT*,'Nombre de points pour lesquels on advect plus que le'
765ccc     &       ,'contenu de la maille : ',n0
766c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
767         DO l=1,llm
768            IF(nl(l).gt.0) THEN
769               iju=0
770c   indicage des mailles concernees par le traitement special
771               DO ij=ijb,ije
772                  IF(iadvplus(ij,l).eq.1.and.mod(ij,iip1).ne.0) THEN
773                     iju=iju+1
774                     indu(iju)=ij
775                  ENDIF
776               ENDDO
777               niju=iju
778c              PRINT*,'niju,nl',niju,nl(l)
779
780c  traitement des mailles
781               DO iju=1,niju
782                  ij=indu(iju)
783                  j=(ij-1)/iip1+1
784                  zu_m=u_m(ij,l)
785                  u_mq(ij,l)=0.
786                  IF(zu_m.gt.0.) THEN
787                     ijq=ij
788                     i=ijq-(j-1)*iip1
789c   accumulation pour les mailles completements advectees
790                     do while(zu_m.gt.masse(ijq,l))
791                        u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l)
792                        zu_m=zu_m-masse(ijq,l)
793                        i=mod(i-2+iim,iim)+1
794                        ijq=(j-1)*iip1+i
795                     ENDDO
796c   ajout de la maille non completement advectee
797                     u_mq(ij,l)=u_mq(ij,l)+zu_m*
798     &               (q(ijq,l)+0.5*(1.-zu_m/masse(ijq,l))*dxq(ijq,l))
799                  ELSE
800                     ijq=ij+1
801                     i=ijq-(j-1)*iip1
802c   accumulation pour les mailles completements advectees
803                     do while(-zu_m.gt.masse(ijq,l))
804                        u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l)
805                        zu_m=zu_m+masse(ijq,l)
806                        i=mod(i,iim)+1
807                        ijq=(j-1)*iip1+i
808                     ENDDO
809c   ajout de la maille non completement advectee
810                     u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l)-
811     &               0.5*(1.+zu_m/masse(ijq,l))*dxq(ijq,l))
812                  ENDIF
813               ENDDO
814            ENDIF
815         ENDDO
816c$OMP END DO NOWAIT
817cym      ENDIF  ! n0.gt.0
818
819
820
821c   bouclage en latitude
822c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
823      DO l=1,llm
824        DO ij=ijb+iip1-1,ije,iip1
825           u_mq(ij,l)=u_mq(ij-iim,l)
826        ENDDO
827      ENDDO
828c$OMP END DO NOWAIT
829
830c   calcul des tendances
831c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
832      DO l=1,llm
833         DO ij=ijb+1,ije
834            new_m=masse(ij,l)+u_m(ij-1,l)-u_m(ij,l)
835            q(ij,l)=(q(ij,l)*masse(ij,l)+
836     &      u_mq(ij-1,l)-u_mq(ij,l))
837     &      /new_m
838            masse(ij,l)=new_m
839         ENDDO
840c   Modif Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
841         DO ij=ijb+iip1-1,ije,iip1
842            q(ij-iim,l)=q(ij,l)
843            masse(ij-iim,l)=masse(ij,l)
844         ENDDO
845      ENDDO
846c$OMP END DO NOWAIT
847c     CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
848c     CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1)
849
850
851      RETURN
852      END
853      SUBROUTINE vlyqs_p(q,pente_max,masse,masse_adv_v,qsat)
854c
855c     Auteurs:   P.Le Van, F.Hourdin, F.Forget 
856c
857c    ********************************************************************
858c     Shema  d'advection " pseudo amont " .
859c    ********************************************************************
860c     q,masse_adv_v,w sont des arguments d'entree  pour le s-pg ....
861c     qsat             est   un argument de sortie pour le s-pg ....
862c
863c
864c   --------------------------------------------------------------------
865      USE parallel_lmdz
866      IMPLICIT NONE
867c
868!-----------------------------------------------------------------------
869!   INCLUDE 'dimensions.h'
870!
871!   dimensions.h contient les dimensions du modele
872!   ndm est tel que iim=2**ndm
873!-----------------------------------------------------------------------
874
875      INTEGER iim,jjm,llm,ndm
876
877      PARAMETER (iim= 128,jjm=96,llm=64,ndm=1)
878
879!-----------------------------------------------------------------------
880!
881! $Header$
882!
883!
884!  ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
885!                 veillez  n'utiliser que des ! pour les commentaires
886!                 et  bien positionner les & des lignes de continuation
887!                 (les placer en colonne 6 et en colonne 73)
888!
889!
890!-----------------------------------------------------------------------
891!   INCLUDE 'paramet.h'
892
893      INTEGER  iip1,iip2,iip3,jjp1,llmp1,llmp2,llmm1
894      INTEGER  kftd,ip1jm,ip1jmp1,ip1jmi1,ijp1llm
895      INTEGER  ijmllm,mvar
896      INTEGER jcfil,jcfllm
897
898      PARAMETER( iip1= iim+1,iip2=iim+2,iip3=iim+3                       &
899     &    ,jjp1=jjm+1-1/jjm)
900      PARAMETER( llmp1 = llm+1,  llmp2 = llm+2, llmm1 = llm-1 )
901      PARAMETER( kftd  = iim/2 -ndm )
902      PARAMETER( ip1jm  = iip1*jjm,  ip1jmp1= iip1*jjp1 )
903      PARAMETER( ip1jmi1= ip1jm - iip1 )
904      PARAMETER( ijp1llm= ip1jmp1 * llm, ijmllm= ip1jm * llm )
905      PARAMETER( mvar= ip1jmp1*( 2*llm+1) + ijmllm )
906      PARAMETER( jcfil=jjm/2+5, jcfllm=jcfil*llm )
907
908!-----------------------------------------------------------------------
909!
910! $Id: logic.h 1520 2011-05-23 11:37:09Z emillour $
911!
912!
913! NB: keep items of different kinds in seperate common blocs to avoid
914!     "misaligned commons" issues
915!-----------------------------------------------------------------------
916! INCLUDE 'logic.h'
917
918      COMMON/logicl/ purmats,forward,leapf,apphys,                      &
919     &  statcl,conser,apdiss,apdelq,saison,ecripar,fxyhypb,ysinus       &
920     &  ,read_start,ok_guide,ok_strato,tidal,ok_gradsfile               &
921     &  ,ok_limit,ok_etat0,hybrid                                       &
922     &  ,moyzon_mu,moyzon_ch
923
924      COMMON/logici/ iflag_phys,iflag_trac
925     
926      LOGICAL purmats,forward,leapf,apphys,statcl,conser,               &
927     & apdiss,apdelq,saison,ecripar,fxyhypb,ysinus                      &
928     &  ,read_start,ok_guide,ok_strato,tidal,ok_gradsfile               &
929     &  ,ok_limit,ok_etat0
930      logical hybrid ! vertical coordinate is hybrid if true (sigma otherwise)
931                     ! (only used if disvert_type==2)
932      logical moyzon_mu,moyzon_ch ! used for zonal averages in Titan
933
934      integer iflag_phys,iflag_trac
935!$OMP THREADPRIVATE(/logicl/)
936!$OMP THREADPRIVATE(/logici/)
937!-----------------------------------------------------------------------
938!
939! $Id: comvert.h 1654 2012-09-24 15:07:18Z aslmd $
940!
941!-----------------------------------------------------------------------
942!   INCLUDE 'comvert.h'
943
944      COMMON/comvertr/ap(llm+1),bp(llm+1),presnivs(llm),dpres(llm),     &
945     &               pa,preff,nivsigs(llm),nivsig(llm+1),               &
946     &               aps(llm),bps(llm),scaleheight,pseudoalt(llm)
947
948      common/comverti/disvert_type, pressure_exner
949
950      real ap     ! hybrid pressure contribution at interlayers
951      real bp     ! hybrid sigma contribution at interlayer
952      real presnivs ! (reference) pressure at mid-layers
953      real dpres
954      real pa     ! reference pressure (Pa) at which hybrid coordinates
955                  ! become purely pressure
956      real preff  ! reference surface pressure (Pa)
957      real nivsigs
958      real nivsig
959      real aps    ! hybrid pressure contribution at mid-layers
960      real bps    ! hybrid sigma contribution at mid-layers
961      real scaleheight ! atmospheric (reference) scale height (km)
962      real pseudoalt ! pseudo-altitude of model levels (km), based on presnivs(),
963                     ! preff and scaleheight
964
965      integer disvert_type ! type of vertical discretization:
966                           ! 1: Earth (default for planet_type==earth),
967                           !     automatic generation
968                           ! 2: Planets (default for planet_type!=earth),
969                           !     using 'z2sig.def' (or 'esasig.def) file
970
971      logical pressure_exner
972!     compute pressure inside layers using Exner function, else use mean
973!     of pressure values at interfaces
974
975 !-----------------------------------------------------------------------
976!
977! $Id: comconst.h 1437 2010-09-30 08:29:10Z emillour $
978!
979!-----------------------------------------------------------------------
980! INCLUDE comconst.h
981
982      COMMON/comconsti/im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl,          &
983     &                 iflag_top_bound,mode_top_bound
984      COMMON/comconstr/dtvr,daysec,                                     &
985     & pi,dtphys,dtdiss,rad,r,kappa,cotot,unsim,g,omeg                  &
986     & ,dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta        &
987     & ,dissip_pupstart  ,tau_top_bound,                                &
988     & daylen,molmass, ihf
989      COMMON/cpdetvenus/cpp,nu_venus,t0_venus
990
991      INTEGER im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl
992      REAL dtvr ! dynamical time step (in s)
993      REAL daysec !length (in s) of a standard day
994      REAL pi    ! something like 3.14159....
995      REAL dtphys ! (s) time step for the physics
996      REAL dtdiss ! (s) time step for the dissipation
997      REAL rad ! (m) radius of the planet
998      REAL r ! Reduced Gas constant r=R/mu
999             ! with R=8.31.. J.K-1.mol-1, mu: mol mass of atmosphere (kg/mol)
1000      REAL cpp   ! Cp
1001      REAL kappa ! kappa=R/Cp
1002      REAL cotot
1003      REAL unsim ! = 1./iim
1004      REAL g ! (m/s2) gravity
1005      REAL omeg ! (rad/s) rotation rate of the planet
1006! Dissipation factors, for Earth model:
1007      REAL dissip_factz,dissip_zref !dissip_deltaz
1008! Dissipation factors, for other planets:
1009      REAL dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta
1010      REAL dissip_pupstart
1011      INTEGER iflag_top_bound,mode_top_bound
1012      REAL tau_top_bound
1013      REAL daylen ! length of solar day, in 'standard' day length
1014      REAL molmass ! (g/mol) molar mass of the atmosphere
1015
1016      REAL nu_venus,t0_venus ! coeffs needed for Cp(T), Venus atmosphere
1017      REAL ihf  ! (W/m2) intrinsic heat flux for giant planets
1018
1019
1020!-----------------------------------------------------------------------
1021!
1022! $Header$
1023!
1024!CDK comgeom
1025      COMMON/comgeom/                                                   &
1026     & cu(ip1jmp1),cv(ip1jm),unscu2(ip1jmp1),unscv2(ip1jm),             &
1027     & aire(ip1jmp1),airesurg(ip1jmp1),aireu(ip1jmp1),                  &
1028     & airev(ip1jm),unsaire(ip1jmp1),apoln,apols,                       &
1029     & unsairez(ip1jm),airuscv2(ip1jm),airvscu2(ip1jm),                 &
1030     & aireij1(ip1jmp1),aireij2(ip1jmp1),aireij3(ip1jmp1),              &
1031     & aireij4(ip1jmp1),alpha1(ip1jmp1),alpha2(ip1jmp1),                &
1032     & alpha3(ip1jmp1),alpha4(ip1jmp1),alpha1p2(ip1jmp1),               &
1033     & alpha1p4(ip1jmp1),alpha2p3(ip1jmp1),alpha3p4(ip1jmp1),           &
1034     & fext(ip1jm),constang(ip1jmp1),rlatu(jjp1),rlatv(jjm),            &
1035     & rlonu(iip1),rlonv(iip1),cuvsurcv(ip1jm),cvsurcuv(ip1jm),         &
1036     & cvusurcu(ip1jmp1),cusurcvu(ip1jmp1),cuvscvgam1(ip1jm),           &
1037     & cuvscvgam2(ip1jm),cvuscugam1(ip1jmp1),                           &
1038     & cvuscugam2(ip1jmp1),cvscuvgam(ip1jm),cuscvugam(ip1jmp1),         &
1039     & unsapolnga1,unsapolnga2,unsapolsga1,unsapolsga2,                 &
1040     & unsair_gam1(ip1jmp1),unsair_gam2(ip1jmp1),unsairz_gam(ip1jm),    &
1041     & aivscu2gam(ip1jm),aiuscv2gam(ip1jm),xprimu(iip1),xprimv(iip1)
1042
1043!
1044        REAL                                                            &
1045     & cu,cv,unscu2,unscv2,aire,airesurg,aireu,airev,unsaire,apoln     ,&
1046     & apols,unsairez,airuscv2,airvscu2,aireij1,aireij2,aireij3,aireij4,&
1047     & alpha1,alpha2,alpha3,alpha4,alpha1p2,alpha1p4,alpha2p3,alpha3p4 ,&
1048     & fext,constang,rlatu,rlatv,rlonu,rlonv,cuvscvgam1,cuvscvgam2     ,&
1049     & cvuscugam1,cvuscugam2,cvscuvgam,cuscvugam,unsapolnga1,unsapolnga2&
1050     & ,unsapolsga1,unsapolsga2,unsair_gam1,unsair_gam2,unsairz_gam    ,&
1051     & aivscu2gam ,aiuscv2gam,cuvsurcv,cvsurcuv,cvusurcu,cusurcvu,xprimu&
1052     & , xprimv
1053!
1054c
1055c
1056c   Arguments:
1057c   ----------
1058      REAL masse(ip1jmp1,llm),pente_max
1059      REAL masse_adv_v( ip1jm,llm)
1060      REAL q(ip1jmp1,llm)
1061      REAL qsat(ip1jmp1,llm)
1062c
1063c      Local 
1064c   ---------
1065c
1066      INTEGER i,ij,l
1067c
1068      REAL airej2,airejjm,airescb(iim),airesch(iim)
1069      REAL dyq(ip1jmp1,llm),dyqv(ip1jm)
1070      REAL adyqv(ip1jm),dyqmax(ip1jmp1)
1071      REAL qbyv(ip1jm,llm)
1072
1073      REAL qpns,qpsn,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
1074c     REAL newq,oldmasse
1075      Logical first
1076      SAVE first
1077c$OMP THREADPRIVATE(first)
1078      REAL convpn,convps,convmpn,convmps
1079      REAL sinlon(iip1),sinlondlon(iip1)
1080      REAL coslon(iip1),coslondlon(iip1)
1081      SAVE sinlon,coslon,sinlondlon,coslondlon
1082      SAVE airej2,airejjm
1083c$OMP THREADPRIVATE(sinlon,coslon,sinlondlon,coslondlon)
1084c$OMP THREADPRIVATE(airej2,airejjm)
1085c
1086c
1087      REAL      SSUM
1088
1089      DATA first/.true./
1090      INTEGER ijb,ije
1091
1092      IF(first) THEN
1093         PRINT*,'Shema  Amont nouveau  appele dans  Vanleer   '
1094         first=.false.
1095         do i=2,iip1
1096            coslon(i)=cos(rlonv(i))
1097            sinlon(i)=sin(rlonv(i))
1098            coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi
1099            sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
1100         ENDDO
1101         coslon(1)=coslon(iip1)
1102         coslondlon(1)=coslondlon(iip1)
1103         sinlon(1)=sinlon(iip1)
1104         sinlondlon(1)=sinlondlon(iip1)
1105         airej2 = SSUM( iim, aire(iip2), 1 )
1106         airejjm= SSUM( iim, aire(ip1jm -iim), 1 ) 
1107      ENDIF
1108
1109c
1110
1111c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1112      DO l = 1, llm
1113c
1114c   --------------------------------
1115c      CALCUL EN LATITUDE
1116c   --------------------------------
1117
1118c   On commence par calculer la valeur du traceur moyenne sur le premier cercle
1119c   de latitude autour du pole (qpns pour le pole nord et qpsn pour
1120c    le pole nord) qui sera utilisee pour evaluer les pentes au pole.
1121
1122      if (pole_nord) then
1123        DO i = 1, iim
1124          airescb(i) = aire(i+ iip1) * q(i+ iip1,l)
1125        ENDDO
1126        qpns   = SSUM( iim,  airescb ,1 ) / airej2
1127      endif
1128     
1129      if (pole_sud) then
1130        DO i = 1, iim
1131          airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l)
1132        ENDDO
1133        qpsn   = SSUM( iim,  airesch ,1 ) / airejjm
1134      endif
1135
1136
1137c   calcul des pentes aux points v
1138
1139      ijb=ij_begin-2*iip1
1140      ije=ij_end+iip1
1141      if (pole_nord) ijb=ij_begin
1142      if (pole_sud)  ije=ij_end-iip1
1143     
1144      DO ij=ijb,ije
1145         dyqv(ij)=q(ij,l)-q(ij+iip1,l)
1146         adyqv(ij)=abs(dyqv(ij))
1147      ENDDO
1148
1149
1150c   calcul des pentes aux points scalaires
1151
1152      ijb=ij_begin-iip1
1153      ije=ij_end+iip1
1154      if (pole_nord) ijb=ij_begin+iip1
1155      if (pole_sud)  ije=ij_end-iip1
1156     
1157      DO ij=ijb,ije
1158         dyq(ij,l)=.5*(dyqv(ij-iip1)+dyqv(ij))
1159         dyqmax(ij)=min(adyqv(ij-iip1),adyqv(ij))
1160         dyqmax(ij)=pente_max*dyqmax(ij)
1161      ENDDO
1162     
1163      IF (pole_nord) THEN
1164
1165c   calcul des pentes aux poles
1166        DO ij=1,iip1
1167           dyq(ij,l)=qpns-q(ij+iip1,l)
1168        ENDDO
1169
1170c   filtrage de la derivee       
1171        dyn1=0.
1172        dyn2=0.
1173        DO ij=1,iim
1174          dyn1=dyn1+sinlondlon(ij)*dyq(ij,l)
1175          dyn2=dyn2+coslondlon(ij)*dyq(ij,l)
1176        ENDDO
1177        DO ij=1,iip1
1178          dyq(ij,l)=dyn1*sinlon(ij)+dyn2*coslon(ij)
1179        ENDDO
1180
1181c   calcul des pentes limites aux poles
1182        fn=1.
1183        DO ij=1,iim
1184          IF(pente_max*adyqv(ij).lt.abs(dyq(ij,l))) THEN
1185            fn=min(pente_max*adyqv(ij)/abs(dyq(ij,l)),fn)
1186          ENDIF
1187        ENDDO
1188     
1189        DO ij=1,iip1
1190         dyq(ij,l)=fn*dyq(ij,l)
1191        ENDDO
1192         
1193      ENDIF
1194     
1195      IF (pole_sud) THEN
1196
1197        DO ij=1,iip1
1198           dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l)-qpsn
1199        ENDDO
1200
1201        dys1=0.
1202        dys2=0.
1203
1204        DO ij=1,iim
1205          dys1=dys1+sinlondlon(ij)*dyq(ip1jm+ij,l)
1206          dys2=dys2+coslondlon(ij)*dyq(ip1jm+ij,l)
1207        ENDDO
1208
1209        DO ij=1,iip1
1210          dyq(ip1jm+ij,l)=dys1*sinlon(ij)+dys2*coslon(ij)
1211        ENDDO
1212       
1213c   calcul des pentes limites aux poles 
1214        fs=1.
1215        DO ij=1,iim
1216        IF(pente_max*adyqv(ij+ip1jm-iip1).lt.abs(dyq(ij+ip1jm,l))) THEN
1217         fs=min(pente_max*adyqv(ij+ip1jm-iip1)/abs(dyq(ij+ip1jm,l)),fs)
1218        ENDIF
1219        ENDDO
1220   
1221        DO ij=1,iip1
1222         dyq(ip1jm+ij,l)=fs*dyq(ip1jm+ij,l)
1223        ENDDO
1224       
1225      ENDIF
1226
1227
1228CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1229C  En memoire de dIFferents tests sur la 
1230C  limitation des pentes aux poles.
1231CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1232C     PRINT*,dyq(1)
1233C     PRINT*,dyqv(iip1+1)
1234C     appn=abs(dyq(1)/dyqv(iip1+1))
1235C     PRINT*,dyq(ip1jm+1)
1236C     PRINT*,dyqv(ip1jm-iip1+1)
1237C     apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
1238C     DO ij=2,iim
1239C        appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
1240C        apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
1241C     ENDDO
1242C     appn=min(pente_max/appn,1.)
1243C     apps=min(pente_max/apps,1.)
1244C
1245C
1246C   cas ou on a un extremum au pole
1247C
1248C     IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
1249C    &   appn=0.
1250C     IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
1251C    &   dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
1252C    &   apps=0.
1253C
1254C   limitation des pentes aux poles
1255C     DO ij=1,iip1
1256C        dyq(ij)=appn*dyq(ij)
1257C        dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
1258C     ENDDO
1259C
1260C   test
1261C      DO ij=1,iip1
1262C         dyq(iip1+ij)=0.
1263C         dyq(ip1jm+ij-iip1)=0.
1264C      ENDDO
1265C      DO ij=1,ip1jmp1
1266C         dyq(ij)=dyq(ij)*cos(rlatu((ij-1)/iip1+1))
1267C      ENDDO
1268C
1269C changement 10 07 96
1270C     IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
1271C    &   THEN
1272C        DO ij=1,iip1
1273C           dyqmax(ij)=0.
1274C        ENDDO
1275C     ELSE
1276C        DO ij=1,iip1
1277C           dyqmax(ij)=pente_max*abs(dyqv(ij))
1278C        ENDDO
1279C     ENDIF
1280C
1281C     IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
1282C    & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
1283C    &THEN
1284C        DO ij=ip1jm+1,ip1jmp1
1285C           dyqmax(ij)=0.
1286C        ENDDO
1287C     ELSE
1288C        DO ij=ip1jm+1,ip1jmp1
1289C           dyqmax(ij)=pente_max*abs(dyqv(ij-iip1))
1290C        ENDDO
1291C     ENDIF
1292C   fin changement 10 07 96
1293CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1294
1295c   calcul des pentes limitees
1296      ijb=ij_begin-iip1
1297      ije=ij_end+iip1
1298      if (pole_nord) ijb=ij_begin+iip1
1299      if (pole_sud)  ije=ij_end-iip1
1300
1301      DO ij=ijb,ije
1302         IF(dyqv(ij)*dyqv(ij-iip1).gt.0.) THEN
1303            dyq(ij,l)=sign(min(abs(dyq(ij,l)),dyqmax(ij)),dyq(ij,l))
1304         ELSE
1305            dyq(ij,l)=0.
1306         ENDIF
1307      ENDDO
1308
1309      ENDDO
1310c$OMP END DO NOWAIT
1311
1312      ijb=ij_begin-iip1
1313      ije=ij_end
1314      if (pole_nord) ijb=ij_begin
1315      if (pole_sud)  ije=ij_end-iip1
1316
1317c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1318      DO l=1,llm
1319       DO ij=ijb,ije
1320         IF( masse_adv_v(ij,l).GT.0. ) THEN
1321           qbyv(ij,l)= MIN( qsat(ij+iip1,l), q(ij+iip1,l )  +
1322     ,      dyq(ij+iip1,l)*0.5*(1.-masse_adv_v(ij,l)/masse(ij+iip1,l)))
1323         ELSE
1324              qbyv(ij,l)= MIN( qsat(ij,l), q(ij,l) - dyq(ij,l) *
1325     ,                   0.5*(1.+masse_adv_v(ij,l)/masse(ij,l)) )
1326         ENDIF
1327          qbyv(ij,l) = masse_adv_v(ij,l)*qbyv(ij,l)
1328       ENDDO
1329      ENDDO
1330c$OMP END DO NOWAIT
1331
1332      ijb=ij_begin
1333      ije=ij_end
1334      if (pole_nord) ijb=ij_begin+iip1
1335      if (pole_sud)  ije=ij_end-iip1
1336
1337c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
1338      DO l=1,llm
1339         DO ij=ijb,ije
1340            newmasse=masse(ij,l)
1341     &      +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
1342            q(ij,l)=(q(ij,l)*masse(ij,l)+qbyv(ij,l)-qbyv(ij-iip1,l))
1343     &         /newmasse
1344            masse(ij,l)=newmasse
1345         ENDDO
1346c.-. ancienne version
1347
1348         IF (pole_nord) THEN
1349
1350           convpn=SSUM(iim,qbyv(1,l),1)/apoln
1351           convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln
1352           DO ij = 1,iip1
1353              newmasse=masse(ij,l)+convmpn*aire(ij)
1354              q(ij,l)=(q(ij,l)*masse(ij,l)+convpn*aire(ij))/
1355     &                 newmasse
1356              masse(ij,l)=newmasse
1357           ENDDO
1358         
1359         ENDIF
1360         
1361         IF (pole_sud) THEN
1362         
1363           convps  = -SSUM(iim,qbyv(ip1jm-iim,l),1)/apols
1364           convmps = -SSUM(iim,masse_adv_v(ip1jm-iim,l),1)/apols
1365           DO ij = ip1jm+1,ip1jmp1
1366              newmasse=masse(ij,l)+convmps*aire(ij)
1367              q(ij,l)=(q(ij,l)*masse(ij,l)+convps*aire(ij))/
1368     &                 newmasse
1369              masse(ij,l)=newmasse
1370           ENDDO
1371         
1372         ENDIF
1373c.-. fin ancienne version
1374
1375c._. nouvelle version
1376c        convpn=SSUM(iim,qbyv(1,l),1)
1377c        convmpn=ssum(iim,masse_adv_v(1,l),1)
1378c        oldmasse=ssum(iim,masse(1,l),1)
1379c        newmasse=oldmasse+convmpn
1380c        newq=(q(1,l)*oldmasse+convpn)/newmasse
1381c        newmasse=newmasse/apoln
1382c        DO ij = 1,iip1
1383c           q(ij,l)=newq
1384c           masse(ij,l)=newmasse*aire(ij)
1385c        ENDDO
1386c        convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
1387c        convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
1388c        oldmasse=ssum(iim,masse(ip1jm-iim,l),1)
1389c        newmasse=oldmasse+convmps
1390c        newq=(q(ip1jmp1,l)*oldmasse+convps)/newmasse
1391c        newmasse=newmasse/apols
1392c        DO ij = ip1jm+1,ip1jmp1
1393c           q(ij,l)=newq
1394c           masse(ij,l)=newmasse*aire(ij)
1395c        ENDDO
1396c._. fin nouvelle version
1397      ENDDO
1398c$OMP END DO NOWAIT
1399      RETURN
1400      END
Note: See TracBrowser for help on using the repository browser.