source: LMDZ5/trunk/libf/vlspltqs_loc.F @ 1630

Last change on this file since 1630 was 1630, checked in by Laurent Fairhead, 12 years ago

Importation initiale du répertoire dyn3dmem


Initial import of dyn3dmem directory

File size: 19.2 KB
Line 
1      SUBROUTINE vlxqs_loc(q,pente_max,masse,u_m,qsat,ijb_x,ije_x)
2c
3c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
4c
5c    ********************************************************************
6c     Shema  d'advection " pseudo amont " .
7c    ********************************************************************
8c
9c   --------------------------------------------------------------------
10      USE parallel
11      IMPLICIT NONE
12c
13#include "dimensions.h"
14#include "paramet.h"
15#include "logic.h"
16#include "comvert.h"
17#include "comconst.h"
18c
19c
20c   Arguments:
21c   ----------
22      REAL masse(ijb_u:ije_u,llm),pente_max
23      REAL u_m( ijb_u:ije_u,llm )
24      REAL q(ijb_u:ije_u,llm)
25      REAL qsat(ijb_u:ije_u,llm)
26c
27c      Local
28c   ---------
29c
30      INTEGER ij,l,j,i,iju,ijq,indu(ijnb_u),niju
31      INTEGER n0,iadvplus(ijb_u:ije_u,llm),nl(llm)
32c
33      REAL new_m,zu_m,zdum(ijb_u:ije_u,llm)
34      REAL dxq(ijb_u:ije_u,llm),dxqu(ijb_u:ije_u)
35      REAL zz(ijb_u:ije_u)
36      REAL adxqu(ijb_u:ije_u),dxqmax(ijb_u:ije_u,llm)
37      REAL u_mq(ijb_u:ije_u,llm)
38
39      REAL      SSUM
40
41
42      INTEGER ijb,ije,ijb_x,ije_x
43     
44
45c   calcul de la pente a droite et a gauche de la maille
46
47c      ijb=ij_begin
48c      ije=ij_end
49
50      ijb=ijb_x
51      ije=ije_x
52       
53      if (pole_nord.and.ijb==1) ijb=ijb+iip1
54      if (pole_sud.and.ije==ip1jmp1)  ije=ije-iip1
55     
56      IF (pente_max.gt.-1.e-5) THEN
57c     IF (pente_max.gt.10) THEN
58
59c   calcul des pentes avec limitation, Van Leer scheme I:
60c   -----------------------------------------------------
61
62c   calcul de la pente aux points u
63
64c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
65         DO l = 1, llm
66            DO ij=ijb,ije-1
67               dxqu(ij)=q(ij+1,l)-q(ij,l)
68c              IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
69c              sigu(ij)=u_m(ij,l)/masse(ij,l)
70            ENDDO
71            DO ij=ijb+iip1-1,ije,iip1
72               dxqu(ij)=dxqu(ij-iim)
73c              sigu(ij)=sigu(ij-iim)
74            ENDDO
75
76            DO ij=ijb,ije
77               adxqu(ij)=abs(dxqu(ij))
78            ENDDO
79
80c   calcul de la pente maximum dans la maille en valeur absolue
81
82            DO ij=ijb+1,ije
83               dxqmax(ij,l)=pente_max*
84     ,      min(adxqu(ij-1),adxqu(ij))
85c limitation subtile
86c    ,      min(adxqu(ij-1)/sigu(ij-1),adxqu(ij)/(1.-sigu(ij)))
87         
88
89            ENDDO
90
91            DO ij=ijb+iip1-1,ije,iip1
92               dxqmax(ij-iim,l)=dxqmax(ij,l)
93            ENDDO
94
95            DO ij=ijb+1,ije
96#ifdef CRAY
97               dxq(ij,l)=
98     ,         cvmgp(dxqu(ij-1)+dxqu(ij),0.,dxqu(ij-1)*dxqu(ij))
99#else
100               IF(dxqu(ij-1)*dxqu(ij).gt.0) THEN
101                  dxq(ij,l)=dxqu(ij-1)+dxqu(ij)
102               ELSE
103c   extremum local
104                  dxq(ij,l)=0.
105               ENDIF
106#endif
107               dxq(ij,l)=0.5*dxq(ij,l)
108               dxq(ij,l)=
109     ,         sign(min(abs(dxq(ij,l)),dxqmax(ij,l)),dxq(ij,l))
110            ENDDO
111
112         ENDDO ! l=1,llm
113c$OMP END DO NOWAIT
114
115      ELSE ! (pente_max.lt.-1.e-5)
116
117c   Pentes produits:
118c   ----------------
119c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
120         DO l = 1, llm
121            DO ij=ijb,ije-1
122               dxqu(ij)=q(ij+1,l)-q(ij,l)
123            ENDDO
124            DO ij=ijb+iip1-1,ije,iip1
125               dxqu(ij)=dxqu(ij-iim)
126            ENDDO
127
128            DO ij=ijb+1,ije
129               zz(ij)=dxqu(ij-1)*dxqu(ij)
130               zz(ij)=zz(ij)+zz(ij)
131               IF(zz(ij).gt.0) THEN
132                  dxq(ij,l)=zz(ij)/(dxqu(ij-1)+dxqu(ij))
133               ELSE
134c   extremum local
135                  dxq(ij,l)=0.
136               ENDIF
137            ENDDO
138
139         ENDDO
140c$OMP END DO NOWAIT
141      ENDIF ! (pente_max.lt.-1.e-5)
142
143c   bouclage de la pente en iip1:
144c   -----------------------------
145c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
146      DO l=1,llm
147         DO ij=ijb+iip1-1,ije,iip1
148            dxq(ij-iim,l)=dxq(ij,l)
149         ENDDO
150
151         DO ij=ijb,ije
152            iadvplus(ij,l)=0
153         ENDDO
154
155      ENDDO
156c$OMP END DO NOWAIT
157     
158      if (pole_nord) THEN
159c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
160        DO l=1,llm     
161          iadvplus(1:iip1,l)=0
162        ENDDO
163c$OMP END DO NOWAIT
164      endif
165     
166      if (pole_sud)  THEN
167c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
168        DO l=1,llm
169          iadvplus(ip1jm+1:ip1jmp1,l)=0
170        ENDDO
171c$OMP END DO NOWAIT
172      endif
173       
174c   calcul des flux a gauche et a droite
175
176#ifdef CRAY
177c--pas encore modification sur Qsat
178c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
179      DO l=1,llm
180       DO ij=ijb,ije-1
181          zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l),
182     ,                     1.+u_m(ij,l)/masse(ij+1,l),
183     ,                     u_m(ij,l))
184          zdum(ij,l)=0.5*zdum(ij,l)
185          u_mq(ij,l)=cvmgp(
186     ,                q(ij,l)+zdum(ij,l)*dxq(ij,l),
187     ,                q(ij+1,l)-zdum(ij,l)*dxq(ij+1,l),
188     ,                u_m(ij,l))
189          u_mq(ij,l)=u_m(ij,l)*u_mq(ij,l)
190       ENDDO
191      ENDDO
192c$OMP END DO NOWAIT
193
194#else
195c   on cumule le flux correspondant a toutes les mailles dont la masse
196c   au travers de la paroi pENDant le pas de temps.
197c   le rapport de melange de l'air advecte est min(q_vanleer, Qsat_downwind)
198c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
199      DO l=1,llm
200       DO ij=ijb,ije-1
201          IF (u_m(ij,l).gt.0.) THEN
202             zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l)
203             u_mq(ij,l)=u_m(ij,l)*
204     $         min(q(ij,l)+0.5*zdum(ij,l)*dxq(ij,l),qsat(ij+1,l))
205          ELSE
206             zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l)
207             u_mq(ij,l)=u_m(ij,l)*
208     $         min(q(ij+1,l)-0.5*zdum(ij,l)*dxq(ij+1,l),qsat(ij,l))
209          ENDIF
210       ENDDO
211      ENDDO
212c$OMP END DO NOWAIT
213#endif
214
215
216c   detection des points ou on advecte plus que la masse de la
217c   maille
218c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
219      DO l=1,llm
220         DO ij=ijb,ije-1
221            IF(zdum(ij,l).lt.0) THEN
222               iadvplus(ij,l)=1
223               u_mq(ij,l)=0.
224            ENDIF
225         ENDDO
226      ENDDO
227c$OMP END DO NOWAIT
228
229c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
230      DO l=1,llm
231       DO ij=ijb+iip1-1,ije,iip1
232          iadvplus(ij,l)=iadvplus(ij-iim,l)
233       ENDDO
234      ENDDO
235c$OMP END DO NOWAIT
236
237
238
239c   traitement special pour le cas ou on advecte en longitude plus que le
240c   contenu de la maille.
241c   cette partie est mal vectorisee.
242
243c   pas d'influence de la pression saturante (pour l'instant)
244
245c  calcul du nombre de maille sur lequel on advecte plus que la maille.
246
247      n0=0
248c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
249      DO l=1,llm
250         nl(l)=0
251         DO ij=ijb,ije
252            nl(l)=nl(l)+iadvplus(ij,l)
253         ENDDO
254         n0=n0+nl(l)
255      ENDDO
256c$OMP END DO NOWAIT
257
258cym ATTENTION ICI en OpenMP reduction pas forcement nécessaire
259cym      IF(n0.gt.1) THEN
260cym        IF(n0.gt.0) THEN
261ccc      PRINT*,'Nombre de points pour lesquels on advect plus que le'
262ccc     &       ,'contenu de la maille : ',n0
263c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
264         DO l=1,llm
265            IF(nl(l).gt.0) THEN
266               iju=0
267c   indicage des mailles concernees par le traitement special
268               DO ij=ijb,ije
269                  IF(iadvplus(ij,l).eq.1.and.mod(ij,iip1).ne.0) THEN
270                     iju=iju+1
271                     indu(iju)=ij
272                  ENDIF
273               ENDDO
274               niju=iju
275c              PRINT*,'niju,nl',niju,nl(l)
276
277c  traitement des mailles
278               DO iju=1,niju
279                  ij=indu(iju)
280                  j=(ij-1)/iip1+1
281                  zu_m=u_m(ij,l)
282                  u_mq(ij,l)=0.
283                  IF(zu_m.gt.0.) THEN
284                     ijq=ij
285                     i=ijq-(j-1)*iip1
286c   accumulation pour les mailles completements advectees
287                     do while(zu_m.gt.masse(ijq,l))
288                        u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l)
289                        zu_m=zu_m-masse(ijq,l)
290                        i=mod(i-2+iim,iim)+1
291                        ijq=(j-1)*iip1+i
292                     ENDDO
293c   ajout de la maille non completement advectee
294                     u_mq(ij,l)=u_mq(ij,l)+zu_m*
295     &               (q(ijq,l)+0.5*(1.-zu_m/masse(ijq,l))*dxq(ijq,l))
296                  ELSE
297                     ijq=ij+1
298                     i=ijq-(j-1)*iip1
299c   accumulation pour les mailles completements advectees
300                     do while(-zu_m.gt.masse(ijq,l))
301                        u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l)
302                        zu_m=zu_m+masse(ijq,l)
303                        i=mod(i,iim)+1
304                        ijq=(j-1)*iip1+i
305                     ENDDO
306c   ajout de la maille non completement advectee
307                     u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l)-
308     &               0.5*(1.+zu_m/masse(ijq,l))*dxq(ijq,l))
309                  ENDIF
310               ENDDO
311            ENDIF
312         ENDDO
313c$OMP END DO NOWAIT
314cym      ENDIF  ! n0.gt.0
315
316
317
318c   bouclage en latitude
319c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
320      DO l=1,llm
321        DO ij=ijb+iip1-1,ije,iip1
322           u_mq(ij,l)=u_mq(ij-iim,l)
323        ENDDO
324      ENDDO
325c$OMP END DO NOWAIT
326
327c   calcul des tendances
328c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
329      DO l=1,llm
330         DO ij=ijb+1,ije
331            new_m=masse(ij,l)+u_m(ij-1,l)-u_m(ij,l)
332            q(ij,l)=(q(ij,l)*masse(ij,l)+
333     &      u_mq(ij-1,l)-u_mq(ij,l))
334     &      /new_m
335            masse(ij,l)=new_m
336         ENDDO
337c   Modif Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
338         DO ij=ijb+iip1-1,ije,iip1
339            q(ij-iim,l)=q(ij,l)
340            masse(ij-iim,l)=masse(ij,l)
341         ENDDO
342      ENDDO
343c$OMP END DO NOWAIT
344c     CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
345c     CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1)
346
347
348      RETURN
349      END
350      SUBROUTINE vlyqs_loc(q,pente_max,masse,masse_adv_v,qsat)
351c
352c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
353c
354c    ********************************************************************
355c     Shema  d'advection " pseudo amont " .
356c    ********************************************************************
357c     q,masse_adv_v,w sont des arguments d'entree  pour le s-pg ....
358c     qsat             est   un argument de sortie pour le s-pg ....
359c
360c
361c   --------------------------------------------------------------------
362      USE parallel
363      IMPLICIT NONE
364c
365#include "dimensions.h"
366#include "paramet.h"
367#include "logic.h"
368#include "comvert.h"
369#include "comconst.h"
370#include "comgeom.h"
371c
372c
373c   Arguments:
374c   ----------
375      REAL masse(ijb_u:ije_u,llm),pente_max
376      REAL masse_adv_v( ijb_v:ije_v,llm)
377      REAL q(ijb_u:ije_u,llm)
378      REAL qsat(ijb_u:ije_u,llm)
379c
380c      Local
381c   ---------
382c
383      INTEGER i,ij,l
384c
385      REAL airej2,airejjm,airescb(iim),airesch(iim)
386      REAL dyq(ijb_u:ije_u,llm),dyqv(ijb_v:ije_v)
387      REAL adyqv(ijb_v:ije_v),dyqmax(ijb_u:ije_u)
388      REAL qbyv(ijb_v:ije_v,llm)
389
390      REAL qpns,qpsn,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
391c     REAL newq,oldmasse
392      Logical first
393      SAVE first
394c$OMP THREADPRIVATE(first)
395      REAL convpn,convps,convmpn,convmps
396      REAL sinlon(iip1),sinlondlon(iip1)
397      REAL coslon(iip1),coslondlon(iip1)
398      SAVE sinlon,coslon,sinlondlon,coslondlon
399      SAVE airej2,airejjm
400c$OMP THREADPRIVATE(sinlon,coslon,sinlondlon,coslondlon)
401c$OMP THREADPRIVATE(airej2,airejjm)
402c
403c
404      REAL      SSUM
405
406      DATA first/.true./
407      INTEGER ijb,ije
408
409      IF(first) THEN
410         PRINT*,'Shema  Amont nouveau  appele dans  Vanleer   '
411         first=.false.
412         do i=2,iip1
413            coslon(i)=cos(rlonv(i))
414            sinlon(i)=sin(rlonv(i))
415            coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi
416            sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
417         ENDDO
418         coslon(1)=coslon(iip1)
419         coslondlon(1)=coslondlon(iip1)
420         sinlon(1)=sinlon(iip1)
421         sinlondlon(1)=sinlondlon(iip1)
422         airej2 = SSUM( iim, aire(iip2), 1 )
423         airejjm= SSUM( iim, aire(ip1jm -iim), 1 )
424      ENDIF
425
426c
427
428c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
429      DO l = 1, llm
430c
431c   --------------------------------
432c      CALCUL EN LATITUDE
433c   --------------------------------
434
435c   On commence par calculer la valeur du traceur moyenne sur le premier cercle
436c   de latitude autour du pole (qpns pour le pole nord et qpsn pour
437c    le pole nord) qui sera utilisee pour evaluer les pentes au pole.
438
439      if (pole_nord) then
440        DO i = 1, iim
441          airescb(i) = aire(i+ iip1) * q(i+ iip1,l)
442        ENDDO
443        qpns   = SSUM( iim,  airescb ,1 ) / airej2
444      endif
445     
446      if (pole_sud) then
447        DO i = 1, iim
448          airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l)
449        ENDDO
450        qpsn   = SSUM( iim,  airesch ,1 ) / airejjm
451      endif
452
453
454c   calcul des pentes aux points v
455
456      ijb=ij_begin-2*iip1
457      ije=ij_end+iip1
458      if (pole_nord) ijb=ij_begin
459      if (pole_sud)  ije=ij_end-iip1
460     
461      DO ij=ijb,ije
462         dyqv(ij)=q(ij,l)-q(ij+iip1,l)
463         adyqv(ij)=abs(dyqv(ij))
464      ENDDO
465
466
467c   calcul des pentes aux points scalaires
468
469      ijb=ij_begin-iip1
470      ije=ij_end+iip1
471      if (pole_nord) ijb=ij_begin+iip1
472      if (pole_sud)  ije=ij_end-iip1
473     
474      DO ij=ijb,ije
475         dyq(ij,l)=.5*(dyqv(ij-iip1)+dyqv(ij))
476         dyqmax(ij)=min(adyqv(ij-iip1),adyqv(ij))
477         dyqmax(ij)=pente_max*dyqmax(ij)
478      ENDDO
479     
480      IF (pole_nord) THEN
481
482c   calcul des pentes aux poles
483        DO ij=1,iip1
484           dyq(ij,l)=qpns-q(ij+iip1,l)
485        ENDDO
486
487c   filtrage de la derivee       
488        dyn1=0.
489        dyn2=0.
490        DO ij=1,iim
491          dyn1=dyn1+sinlondlon(ij)*dyq(ij,l)
492          dyn2=dyn2+coslondlon(ij)*dyq(ij,l)
493        ENDDO
494        DO ij=1,iip1
495          dyq(ij,l)=dyn1*sinlon(ij)+dyn2*coslon(ij)
496        ENDDO
497
498c   calcul des pentes limites aux poles
499        fn=1.
500        DO ij=1,iim
501          IF(pente_max*adyqv(ij).lt.abs(dyq(ij,l))) THEN
502            fn=min(pente_max*adyqv(ij)/abs(dyq(ij,l)),fn)
503          ENDIF
504        ENDDO
505     
506        DO ij=1,iip1
507         dyq(ij,l)=fn*dyq(ij,l)
508        ENDDO
509         
510      ENDIF
511     
512      IF (pole_sud) THEN
513
514        DO ij=1,iip1
515           dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l)-qpsn
516        ENDDO
517
518        dys1=0.
519        dys2=0.
520
521        DO ij=1,iim
522          dys1=dys1+sinlondlon(ij)*dyq(ip1jm+ij,l)
523          dys2=dys2+coslondlon(ij)*dyq(ip1jm+ij,l)
524        ENDDO
525
526        DO ij=1,iip1
527          dyq(ip1jm+ij,l)=dys1*sinlon(ij)+dys2*coslon(ij)
528        ENDDO
529       
530c   calcul des pentes limites aux poles
531        fs=1.
532        DO ij=1,iim
533        IF(pente_max*adyqv(ij+ip1jm-iip1).lt.abs(dyq(ij+ip1jm,l))) THEN
534         fs=min(pente_max*adyqv(ij+ip1jm-iip1)/abs(dyq(ij+ip1jm,l)),fs)
535        ENDIF
536        ENDDO
537   
538        DO ij=1,iip1
539         dyq(ip1jm+ij,l)=fs*dyq(ip1jm+ij,l)
540        ENDDO
541       
542      ENDIF
543
544
545CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
546C  En memoire de dIFferents tests sur la
547C  limitation des pentes aux poles.
548CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
549C     PRINT*,dyq(1)
550C     PRINT*,dyqv(iip1+1)
551C     apn=abs(dyq(1)/dyqv(iip1+1))
552C     PRINT*,dyq(ip1jm+1)
553C     PRINT*,dyqv(ip1jm-iip1+1)
554C     aps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
555C     DO ij=2,iim
556C        apn=amax1(abs(dyq(ij)/dyqv(ij)),apn)
557C        aps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),aps)
558C     ENDDO
559C     apn=min(pente_max/apn,1.)
560C     aps=min(pente_max/aps,1.)
561C
562C
563C   cas ou on a un extremum au pole
564C
565C     IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
566C    &   apn=0.
567C     IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
568C    &   dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
569C    &   aps=0.
570C
571C   limitation des pentes aux poles
572C     DO ij=1,iip1
573C        dyq(ij)=apn*dyq(ij)
574C        dyq(ip1jm+ij)=aps*dyq(ip1jm+ij)
575C     ENDDO
576C
577C   test
578C      DO ij=1,iip1
579C         dyq(iip1+ij)=0.
580C         dyq(ip1jm+ij-iip1)=0.
581C      ENDDO
582C      DO ij=1,ip1jmp1
583C         dyq(ij)=dyq(ij)*cos(rlatu((ij-1)/iip1+1))
584C      ENDDO
585C
586C changement 10 07 96
587C     IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
588C    &   THEN
589C        DO ij=1,iip1
590C           dyqmax(ij)=0.
591C        ENDDO
592C     ELSE
593C        DO ij=1,iip1
594C           dyqmax(ij)=pente_max*abs(dyqv(ij))
595C        ENDDO
596C     ENDIF
597C
598C     IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
599C    & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
600C    &THEN
601C        DO ij=ip1jm+1,ip1jmp1
602C           dyqmax(ij)=0.
603C        ENDDO
604C     ELSE
605C        DO ij=ip1jm+1,ip1jmp1
606C           dyqmax(ij)=pente_max*abs(dyqv(ij-iip1))
607C        ENDDO
608C     ENDIF
609C   fin changement 10 07 96
610CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
611
612c   calcul des pentes limitees
613      ijb=ij_begin-iip1
614      ije=ij_end+iip1
615      if (pole_nord) ijb=ij_begin+iip1
616      if (pole_sud)  ije=ij_end-iip1
617
618      DO ij=ijb,ije
619         IF(dyqv(ij)*dyqv(ij-iip1).gt.0.) THEN
620            dyq(ij,l)=sign(min(abs(dyq(ij,l)),dyqmax(ij)),dyq(ij,l))
621         ELSE
622            dyq(ij,l)=0.
623         ENDIF
624      ENDDO
625
626      ENDDO
627c$OMP END DO NOWAIT
628
629      ijb=ij_begin-iip1
630      ije=ij_end
631      if (pole_nord) ijb=ij_begin
632      if (pole_sud)  ije=ij_end-iip1
633
634c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
635      DO l=1,llm
636       DO ij=ijb,ije
637         IF( masse_adv_v(ij,l).GT.0. ) THEN
638           qbyv(ij,l)= MIN( qsat(ij+iip1,l), q(ij+iip1,l )  +
639     ,      dyq(ij+iip1,l)*0.5*(1.-masse_adv_v(ij,l)/masse(ij+iip1,l)))
640         ELSE
641              qbyv(ij,l)= MIN( qsat(ij,l), q(ij,l) - dyq(ij,l) *
642     ,                   0.5*(1.+masse_adv_v(ij,l)/masse(ij,l)) )
643         ENDIF
644          qbyv(ij,l) = masse_adv_v(ij,l)*qbyv(ij,l)
645       ENDDO
646      ENDDO
647c$OMP END DO NOWAIT
648
649      ijb=ij_begin
650      ije=ij_end
651      if (pole_nord) ijb=ij_begin+iip1
652      if (pole_sud)  ije=ij_end-iip1
653
654c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
655      DO l=1,llm
656         DO ij=ijb,ije
657            newmasse=masse(ij,l)
658     &      +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
659            q(ij,l)=(q(ij,l)*masse(ij,l)+qbyv(ij,l)-qbyv(ij-iip1,l))
660     &         /newmasse
661            masse(ij,l)=newmasse
662         ENDDO
663c.-. ancienne version
664
665         IF (pole_nord) THEN
666
667           convpn=SSUM(iim,qbyv(1,l),1)/apoln
668           convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln
669           DO ij = 1,iip1
670              newmasse=masse(ij,l)+convmpn*aire(ij)
671              q(ij,l)=(q(ij,l)*masse(ij,l)+convpn*aire(ij))/
672     &                 newmasse
673              masse(ij,l)=newmasse
674           ENDDO
675         
676         ENDIF
677         
678         IF (pole_sud) THEN
679         
680           convps  = -SSUM(iim,qbyv(ip1jm-iim,l),1)/apols
681           convmps = -SSUM(iim,masse_adv_v(ip1jm-iim,l),1)/apols
682           DO ij = ip1jm+1,ip1jmp1
683              newmasse=masse(ij,l)+convmps*aire(ij)
684              q(ij,l)=(q(ij,l)*masse(ij,l)+convps*aire(ij))/
685     &                 newmasse
686              masse(ij,l)=newmasse
687           ENDDO
688         
689         ENDIF
690c.-. fin ancienne version
691
692c._. nouvelle version
693c        convpn=SSUM(iim,qbyv(1,l),1)
694c        convmpn=ssum(iim,masse_adv_v(1,l),1)
695c        oldmasse=ssum(iim,masse(1,l),1)
696c        newmasse=oldmasse+convmpn
697c        newq=(q(1,l)*oldmasse+convpn)/newmasse
698c        newmasse=newmasse/apoln
699c        DO ij = 1,iip1
700c           q(ij,l)=newq
701c           masse(ij,l)=newmasse*aire(ij)
702c        ENDDO
703c        convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
704c        convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
705c        oldmasse=ssum(iim,masse(ip1jm-iim,l),1)
706c        newmasse=oldmasse+convmps
707c        newq=(q(ip1jmp1,l)*oldmasse+convps)/newmasse
708c        newmasse=newmasse/apols
709c        DO ij = ip1jm+1,ip1jmp1
710c           q(ij,l)=newq
711c           masse(ij,l)=newmasse*aire(ij)
712c        ENDDO
713c._. fin nouvelle version
714      ENDDO
715c$OMP END DO NOWAIT
716      RETURN
717      END
Note: See TracBrowser for help on using the repository browser.