source: LMDZ5/trunk/libf/vlsplt_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: 24.1 KB
Line 
1      SUBROUTINE vlx_loc(q,pente_max,masse,u_m,ijb_x,ije_x)
2
3c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
4c
5c    ********************************************************************
6c     Shema  d'advection " pseudo amont " .
7c    ********************************************************************
8c     nq,iq,q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....
9c
10c
11c   --------------------------------------------------------------------
12      USE Parallel
13      IMPLICIT NONE
14c
15#include "dimensions.h"
16#include "paramet.h"
17#include "logic.h"
18#include "comvert.h"
19#include "comconst.h"
20c
21c
22c   Arguments:
23c   ----------
24      REAL masse(ijb_u:ije_u,llm),pente_max
25      REAL u_m( ijb_u:ije_u,llm ),pbarv( iip1,jjb_v:jje_v,llm)
26      REAL q(ijb_u:ije_u,llm)
27      REAL w(ijb_u:ije_u,llm)
28c
29c      Local
30c   ---------
31c
32      INTEGER ij,l,j,i,iju,ijq,indu(ijnb_u),niju
33      INTEGER n0,iadvplus(ijb_u:ije_u,llm),nl(llm)
34c
35      REAL new_m,zu_m,zdum(ijb_u:ije_u,llm)
36      REAL sigu(ijb_u:ije_u),dxq(ijb_u:ije_u,llm),dxqu(ijb_u:ije_u)
37      REAL zz(ijb_u:ije_u)
38      REAL adxqu(ijb_u:ije_u),dxqmax(ijb_u:ije_u,llm)
39      REAL u_mq(ijb_u:ije_u,llm)
40
41      Logical extremum
42
43      REAL      SSUM
44      EXTERNAL  SSUM
45
46      REAL z1,z2,z3
47
48      INTEGER ijb,ije,ijb_x,ije_x
49     
50c   calcul de la pente a droite et a gauche de la maille
51
52      ijb=ijb_x
53      ije=ije_x
54       
55      if (pole_nord.and.ijb==1) ijb=ijb+iip1
56      if (pole_sud.and.ije==ip1jmp1)  ije=ije-iip1
57         
58      IF (pente_max.gt.-1.e-5) THEN
59c       IF (pente_max.gt.10) THEN
60
61c   calcul des pentes avec limitation, Van Leer scheme I:
62c   -----------------------------------------------------
63
64c   calcul de la pente aux points u
65c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)         
66         DO l = 1, llm
67           
68            DO ij=ijb,ije-1
69               dxqu(ij)=q(ij+1,l)-q(ij,l)
70c              IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
71c              sigu(ij)=u_m(ij,l)/masse(ij,l)
72            ENDDO
73            DO ij=ijb+iip1-1,ije,iip1
74               dxqu(ij)=dxqu(ij-iim)
75c              sigu(ij)=sigu(ij-iim)
76            ENDDO
77
78            DO ij=ijb,ije
79               adxqu(ij)=abs(dxqu(ij))
80            ENDDO
81
82c   calcul de la pente maximum dans la maille en valeur absolue
83
84            DO ij=ijb+1,ije
85               dxqmax(ij,l)=pente_max*
86     ,      min(adxqu(ij-1),adxqu(ij))
87c limitation subtile
88c    ,      min(adxqu(ij-1)/sigu(ij-1),adxqu(ij)/(1.-sigu(ij)))
89         
90
91            ENDDO
92
93            DO ij=ijb+iip1-1,ije,iip1
94               dxqmax(ij-iim,l)=dxqmax(ij,l)
95            ENDDO
96
97            DO ij=ijb+1,ije
98#ifdef CRAY
99               dxq(ij,l)=
100     ,         cvmgp(dxqu(ij-1)+dxqu(ij),0.,dxqu(ij-1)*dxqu(ij))
101#else
102               IF(dxqu(ij-1)*dxqu(ij).gt.0) THEN
103                  dxq(ij,l)=dxqu(ij-1)+dxqu(ij)
104               ELSE
105c   extremum local
106                  dxq(ij,l)=0.
107               ENDIF
108#endif
109               dxq(ij,l)=0.5*dxq(ij,l)
110               dxq(ij,l)=
111     ,         sign(min(abs(dxq(ij,l)),dxqmax(ij,l)),dxq(ij,l))
112            ENDDO
113
114         ENDDO ! l=1,llm
115c$OMP END DO NOWAIT
116c       print*,'Ok calcul des pentes'
117
118      ELSE ! (pente_max.lt.-1.e-5)
119
120c   Pentes produits:
121c   ----------------
122c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
123         DO l = 1, llm
124            DO ij=ijb,ije-1
125               dxqu(ij)=q(ij+1,l)-q(ij,l)
126            ENDDO
127            DO ij=ijb+iip1-1,ije,iip1
128               dxqu(ij)=dxqu(ij-iim)
129            ENDDO
130
131            DO ij=ijb+1,ije
132               zz(ij)=dxqu(ij-1)*dxqu(ij)
133               zz(ij)=zz(ij)+zz(ij)
134               IF(zz(ij).gt.0) THEN
135                  dxq(ij,l)=zz(ij)/(dxqu(ij-1)+dxqu(ij))
136               ELSE
137c   extremum local
138                  dxq(ij,l)=0.
139               ENDIF
140            ENDDO
141
142         ENDDO
143c$OMP END DO NOWAIT
144      ENDIF ! (pente_max.lt.-1.e-5)
145
146c   bouclage de la pente en iip1:
147c   -----------------------------
148c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
149      DO l=1,llm
150         DO ij=ijb+iip1-1,ije,iip1
151            dxq(ij-iim,l)=dxq(ij,l)
152         ENDDO
153         DO ij=ijb,ije
154            iadvplus(ij,l)=0
155         ENDDO
156
157      ENDDO
158c$OMP END DO NOWAIT
159c        print*,'Bouclage en iip1'
160
161c   calcul des flux a gauche et a droite
162
163#ifdef CRAY
164c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
165      DO l=1,llm
166       DO ij=ijb,ije-1
167          zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l),
168     ,                     1.+u_m(ij,l)/masse(ij+1,l),
169     ,                     u_m(ij,l))
170          zdum(ij,l)=0.5*zdum(ij,l)
171          u_mq(ij,l)=cvmgp(
172     ,                q(ij,l)+zdum(ij,l)*dxq(ij,l),
173     ,                q(ij+1,l)-zdum(ij,l)*dxq(ij+1,l),
174     ,                u_m(ij,l))
175          u_mq(ij,l)=u_m(ij,l)*u_mq(ij,l)
176       ENDDO
177      ENDDO
178c$OMP END DO NOWAIT
179#else
180c   on cumule le flux correspondant a toutes les mailles dont la masse
181c   au travers de la paroi pENDant le pas de temps.
182c       print*,'Cumule ....'
183c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
184      DO l=1,llm
185       DO ij=ijb,ije-1
186c       print*,'masse(',ij,')=',masse(ij,l)
187          IF (u_m(ij,l).gt.0.) THEN
188             zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l)
189             u_mq(ij,l)=u_m(ij,l)*(q(ij,l)+0.5*zdum(ij,l)*dxq(ij,l))
190          ELSE
191             zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l)
192             u_mq(ij,l)=u_m(ij,l)*(q(ij+1,l)-0.5*zdum(ij,l)*dxq(ij+1,l))
193          ENDIF
194       ENDDO
195      ENDDO
196c$OMP END DO NOWAIT
197#endif
198c       stop
199
200c       go to 9999
201c   detection des points ou on advecte plus que la masse de la
202c   maille
203c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
204      DO l=1,llm
205         DO ij=ijb,ije-1
206            IF(zdum(ij,l).lt.0) THEN
207               iadvplus(ij,l)=1
208               u_mq(ij,l)=0.
209            ENDIF
210         ENDDO
211      ENDDO
212c$OMP END DO NOWAIT
213c       print*,'Ok test 1'
214c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
215      DO l=1,llm
216       DO ij=ijb+iip1-1,ije,iip1
217          iadvplus(ij,l)=iadvplus(ij-iim,l)
218       ENDDO
219      ENDDO
220c$OMP END DO NOWAIT
221c        print*,'Ok test 2'
222
223
224c   traitement special pour le cas ou on advecte en longitude plus que le
225c   contenu de la maille.
226c   cette partie est mal vectorisee.
227
228c  calcul du nombre de maille sur lequel on advecte plus que la maille.
229
230      n0=0
231c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
232      DO l=1,llm
233         nl(l)=0
234         DO ij=ijb,ije
235            nl(l)=nl(l)+iadvplus(ij,l)
236         ENDDO
237         n0=n0+nl(l)
238      ENDDO
239c$OMP END DO NOWAIT
240cym      IF(n0.gt.1) THEN
241cym      IF(n0.gt.0) THEN
242
243c      PRINT*,'Nombre de points pour lesquels on advect plus que le'
244c     &       ,'contenu de la maille : ',n0
245c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
246         DO l=1,llm
247            IF(nl(l).gt.0) THEN
248               iju=0
249c   indicage des mailles concernees par le traitement special
250               DO ij=ijb,ije
251                  IF(iadvplus(ij,l).eq.1.and.mod(ij,iip1).ne.0) THEN
252                     iju=iju+1
253                     indu(iju)=ij
254                  ENDIF
255               ENDDO
256               niju=iju
257c              PRINT*,'niju,nl',niju,nl(l)
258
259c  traitement des mailles
260               DO iju=1,niju
261                  ij=indu(iju)
262                  j=(ij-1)/iip1+1
263                  zu_m=u_m(ij,l)
264                  u_mq(ij,l)=0.
265                  IF(zu_m.gt.0.) THEN
266                     ijq=ij
267                     i=ijq-(j-1)*iip1
268c   accumulation pour les mailles completements advectees
269                     do while(zu_m.gt.masse(ijq,l))
270                        u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l)
271                        zu_m=zu_m-masse(ijq,l)
272                        i=mod(i-2+iim,iim)+1
273                        ijq=(j-1)*iip1+i
274                     ENDDO
275c   ajout de la maille non completement advectee
276                     u_mq(ij,l)=u_mq(ij,l)+zu_m*
277     &               (q(ijq,l)+0.5*(1.-zu_m/masse(ijq,l))*dxq(ijq,l))
278                  ELSE
279                     ijq=ij+1
280                     i=ijq-(j-1)*iip1
281c   accumulation pour les mailles completements advectees
282                     do while(-zu_m.gt.masse(ijq,l))
283                        u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l)
284                        zu_m=zu_m+masse(ijq,l)
285                        i=mod(i,iim)+1
286                        ijq=(j-1)*iip1+i
287                     ENDDO
288c   ajout de la maille non completement advectee
289                     u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l)-
290     &               0.5*(1.+zu_m/masse(ijq,l))*dxq(ijq,l))
291                  ENDIF
292               ENDDO
293            ENDIF
294         ENDDO
295c$OMP END DO NOWAIT
296cym      ENDIF  ! n0.gt.0
2979999    continue
298
299
300c   bouclage en latitude
301c       print*,'Avant bouclage en latitude'
302c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
303      DO l=1,llm
304        DO ij=ijb+iip1-1,ije,iip1
305           u_mq(ij,l)=u_mq(ij-iim,l)
306        ENDDO
307      ENDDO
308c$OMP END DO NOWAIT
309
310c   calcul des tENDances
311c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
312      DO l=1,llm
313         DO ij=ijb+1,ije
314            new_m=masse(ij,l)+u_m(ij-1,l)-u_m(ij,l)
315            q(ij,l)=(q(ij,l)*masse(ij,l)+
316     &      u_mq(ij-1,l)-u_mq(ij,l))
317     &      /new_m
318            masse(ij,l)=new_m
319         ENDDO
320c   ModIF Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
321         DO ij=ijb+iip1-1,ije,iip1
322            q(ij-iim,l)=q(ij,l)
323            masse(ij-iim,l)=masse(ij,l)
324         ENDDO
325      ENDDO
326c$OMP END DO NOWAIT
327c     CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
328c     CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1)
329
330
331      RETURN
332      END
333
334
335      SUBROUTINE vly_loc(q,pente_max,masse,masse_adv_v)
336c
337c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
338c
339c    ********************************************************************
340c     Shema  d'advection " pseudo amont " .
341c    ********************************************************************
342c     q,masse_adv_v,w sont des arguments d'entree  pour le s-pg ....
343c     dq               sont des arguments de sortie pour le s-pg ....
344c
345c
346c   --------------------------------------------------------------------
347      USE parallel
348      IMPLICIT NONE
349c
350#include "dimensions.h"
351#include "paramet.h"
352#include "logic.h"
353#include "comvert.h"
354#include "comconst.h"
355#include "comgeom.h"
356c
357c
358c   Arguments:
359c   ----------
360      REAL masse(ijb_u:ije_u,llm),pente_max
361      REAL masse_adv_v( ijb_v:ije_v,llm)
362      REAL q(ijb_u:ije_u,llm), dq( ijb_u:ije_u,llm)
363c
364c      Local
365c   ---------
366c
367      INTEGER i,ij,l
368c
369      REAL airej2,airejjm,airescb(iim),airesch(iim)
370      REAL dyq(ijb_u:ije_u,llm),dyqv(ijb_v:ije_v),zdvm(ijb_u:ije_u,llm)
371      REAL adyqv(ijb_v:ije_v),dyqmax(ijb_u:ije_u)
372      REAL qbyv(ijb_v:ije_v,llm)
373
374      REAL qpns,qpsn,apn,aps,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
375c     REAL newq,oldmasse
376      Logical extremum,first,testcpu
377      REAL temps0,temps1,temps2,temps3,temps4,temps5,second
378      SAVE temps0,temps1,temps2,temps3,temps4,temps5
379c$OMP THREADPRIVATE(temps0,temps1,temps2,temps3,temps4,temps5)
380      SAVE first,testcpu
381c$OMP THREADPRIVATE(first,testcpu)
382
383      REAL convpn,convps,convmpn,convmps
384      real massepn,masseps,qpn,qps
385      REAL sinlon(iip1),sinlondlon(iip1)
386      REAL coslon(iip1),coslondlon(iip1)
387      SAVE sinlon,coslon,sinlondlon,coslondlon
388c$OMP THREADPRIVATE(sinlon,coslon,sinlondlon,coslondlon)
389      SAVE airej2,airejjm
390c$OMP THREADPRIVATE(airej2,airejjm)
391c
392c
393      REAL      SSUM
394      EXTERNAL  SSUM
395
396      DATA first,testcpu/.true.,.false./
397      DATA temps0,temps1,temps2,temps3,temps4,temps5/0.,0.,0.,0.,0.,0./
398      INTEGER ijb,ije
399
400      IF(first) THEN
401c         PRINT*,'Shema  Amont nouveau  appele dans  Vanleer   '
402         first=.false.
403         do i=2,iip1
404            coslon(i)=cos(rlonv(i))
405            sinlon(i)=sin(rlonv(i))
406            coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi
407            sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
408         ENDDO
409         coslon(1)=coslon(iip1)
410         coslondlon(1)=coslondlon(iip1)
411         sinlon(1)=sinlon(iip1)
412         sinlondlon(1)=sinlondlon(iip1)
413         airej2 = SSUM( iim, aire(iip2), 1 )
414         airejjm= SSUM( iim, aire(ip1jm -iim), 1 )
415      ENDIF
416
417c
418c       PRINT*,'CALCUL EN LATITUDE'
419
420c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
421      DO l = 1, llm
422c
423c   --------------------------------
424c      CALCUL EN LATITUDE
425c   --------------------------------
426
427c   On commence par calculer la valeur du traceur moyenne sur le premier cercle
428c   de latitude autour du pole (qpns pour le pole nord et qpsn pour
429c    le pole nord) qui sera utilisee pour evaluer les pentes au pole.
430     
431      if (pole_nord) then
432        DO i = 1, iim
433          airescb(i) = aire(i+ iip1) * q(i+ iip1,l)
434        ENDDO
435        qpns   = SSUM( iim,  airescb ,1 ) / airej2
436      endif
437     
438      if (pole_sud) then
439        DO i = 1, iim
440          airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l)
441        ENDDO
442        qpsn   = SSUM( iim,  airesch ,1 ) / airejjm
443      endif
444     
445     
446
447c   calcul des pentes aux points v
448
449      ijb=ij_begin-2*iip1
450      ije=ij_end+iip1
451      if (pole_nord) ijb=ij_begin
452      if (pole_sud)  ije=ij_end-iip1
453     
454      DO ij=ijb,ije
455         dyqv(ij)=q(ij,l)-q(ij+iip1,l)
456         adyqv(ij)=abs(dyqv(ij))
457      ENDDO
458
459c   calcul des pentes aux points scalaires
460      ijb=ij_begin-iip1
461      ije=ij_end+iip1
462      if (pole_nord) ijb=ij_begin+iip1
463      if (pole_sud)  ije=ij_end-iip1
464     
465      DO ij=ijb,ije
466         dyq(ij,l)=.5*(dyqv(ij-iip1)+dyqv(ij))
467         dyqmax(ij)=min(adyqv(ij-iip1),adyqv(ij))
468         dyqmax(ij)=pente_max*dyqmax(ij)
469      ENDDO
470
471c   calcul des pentes aux poles
472      IF (pole_nord) THEN
473        DO ij=1,iip1
474           dyq(ij,l)=qpns-q(ij+iip1,l)
475        ENDDO
476       
477        dyn1=0.
478        dyn2=0.
479        DO ij=1,iim
480          dyn1=dyn1+sinlondlon(ij)*dyq(ij,l)
481          dyn2=dyn2+coslondlon(ij)*dyq(ij,l)
482        ENDDO
483        DO ij=1,iip1
484          dyq(ij,l)=dyn1*sinlon(ij)+dyn2*coslon(ij)
485        ENDDO
486       
487        DO ij=1,iip1
488         dyq(ij,l)=0.
489        ENDDO
490c ym tout cela ne sert pas a grand chose
491      ENDIF
492     
493      IF (pole_sud) THEN
494
495        DO ij=1,iip1
496           dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l)-qpsn
497        ENDDO
498
499        dys1=0.
500        dys2=0.
501
502        DO ij=1,iim
503          dys1=dys1+sinlondlon(ij)*dyq(ip1jm+ij,l)
504          dys2=dys2+coslondlon(ij)*dyq(ip1jm+ij,l)
505        ENDDO
506
507        DO ij=1,iip1
508          dyq(ip1jm+ij,l)=dys1*sinlon(ij)+dys2*coslon(ij)
509        ENDDO
510       
511        DO ij=1,iip1
512         dyq(ip1jm+ij,l)=0.
513        ENDDO
514c ym tout cela ne sert pas a grand chose
515      ENDIF
516
517c   filtrage de la derivee
518
519c   calcul des pentes limites aux poles
520c ym partie inutile
521c      goto 8888
522c      fn=1.
523c      fs=1.
524c      DO ij=1,iim
525c         IF(pente_max*adyqv(ij).lt.abs(dyq(ij,l))) THEN
526c            fn=min(pente_max*adyqv(ij)/abs(dyq(ij,l)),fn)
527c         ENDIF
528c      IF(pente_max*adyqv(ij+ip1jm-iip1).lt.abs(dyq(ij+ip1jm,l))) THEN
529c         fs=min(pente_max*adyqv(ij+ip1jm-iip1)/abs(dyq(ij+ip1jm,l)),fs)
530c         ENDIF
531c      ENDDO
532c      DO ij=1,iip1
533c         dyq(ij,l)=fn*dyq(ij,l)
534c         dyq(ip1jm+ij,l)=fs*dyq(ip1jm+ij,l)
535c      ENDDO
536c 8888    continue
537
538
539CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
540C  En memoire de dIFferents tests sur la
541C  limitation des pentes aux poles.
542CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
543C     PRINT*,dyq(1)
544C     PRINT*,dyqv(iip1+1)
545C     apn=abs(dyq(1)/dyqv(iip1+1))
546C     PRINT*,dyq(ip1jm+1)
547C     PRINT*,dyqv(ip1jm-iip1+1)
548C     aps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
549C     DO ij=2,iim
550C        apn=amax1(abs(dyq(ij)/dyqv(ij)),apn)
551C        aps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),aps)
552C     ENDDO
553C     apn=min(pente_max/apn,1.)
554C     aps=min(pente_max/aps,1.)
555C
556C
557C   cas ou on a un extremum au pole
558C
559C     IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
560C    &   apn=0.
561C     IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
562C    &   dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
563C    &   aps=0.
564C
565C   limitation des pentes aux poles
566C     DO ij=1,iip1
567C        dyq(ij)=apn*dyq(ij)
568C        dyq(ip1jm+ij)=aps*dyq(ip1jm+ij)
569C     ENDDO
570C
571C   test
572C      DO ij=1,iip1
573C         dyq(iip1+ij)=0.
574C         dyq(ip1jm+ij-iip1)=0.
575C      ENDDO
576C      DO ij=1,ip1jmp1
577C         dyq(ij)=dyq(ij)*cos(rlatu((ij-1)/iip1+1))
578C      ENDDO
579C
580C changement 10 07 96
581C     IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
582C    &   THEN
583C        DO ij=1,iip1
584C           dyqmax(ij)=0.
585C        ENDDO
586C     ELSE
587C        DO ij=1,iip1
588C           dyqmax(ij)=pente_max*abs(dyqv(ij))
589C        ENDDO
590C     ENDIF
591C
592C     IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
593C    & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
594C    &THEN
595C        DO ij=ip1jm+1,ip1jmp1
596C           dyqmax(ij)=0.
597C        ENDDO
598C     ELSE
599C        DO ij=ip1jm+1,ip1jmp1
600C           dyqmax(ij)=pente_max*abs(dyqv(ij-iip1))
601C        ENDDO
602C     ENDIF
603C   fin changement 10 07 96
604CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
605
606c   calcul des pentes limitees
607      ijb=ij_begin-iip1
608      ije=ij_end+iip1
609      if (pole_nord) ijb=ij_begin+iip1
610      if (pole_sud)  ije=ij_end-iip1
611
612      DO ij=ijb,ije
613         IF(dyqv(ij)*dyqv(ij-iip1).gt.0.) THEN
614            dyq(ij,l)=sign(min(abs(dyq(ij,l)),dyqmax(ij)),dyq(ij,l))
615         ELSE
616            dyq(ij,l)=0.
617         ENDIF
618      ENDDO
619
620      ENDDO
621c$OMP END DO NOWAIT
622
623      ijb=ij_begin-iip1
624      ije=ij_end
625      if (pole_nord) ijb=ij_begin
626      if (pole_sud)  ije=ij_end-iip1
627
628c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
629      DO l=1,llm
630       DO ij=ijb,ije
631          IF(masse_adv_v(ij,l).gt.0) THEN
632              qbyv(ij,l)=q(ij+iip1,l)+dyq(ij+iip1,l)*
633     ,                   0.5*(1.-masse_adv_v(ij,l)/masse(ij+iip1,l))
634          ELSE
635              qbyv(ij,l)=q(ij,l)-dyq(ij,l)*
636     ,                   0.5*(1.+masse_adv_v(ij,l)/masse(ij,l))
637          ENDIF
638          qbyv(ij,l)=masse_adv_v(ij,l)*qbyv(ij,l)
639       ENDDO
640      ENDDO
641c$OMP END DO NOWAIT
642     
643      ijb=ij_begin
644      ije=ij_end
645      if (pole_nord) ijb=ij_begin+iip1
646      if (pole_sud)  ije=ij_end-iip1
647     
648c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
649      DO l=1,llm
650         DO ij=ijb,ije
651            newmasse=masse(ij,l)
652     &      +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
653     
654            q(ij,l)=(q(ij,l)*masse(ij,l)+qbyv(ij,l)-qbyv(ij-iip1,l))
655     &         /newmasse
656            masse(ij,l)=newmasse
657         ENDDO
658c.-. ancienne version
659c        convpn=SSUM(iim,qbyv(1,l),1)/apoln
660c        convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln
661         if (pole_nord) then
662           convpn=SSUM(iim,qbyv(1,l),1)
663           convmpn=ssum(iim,masse_adv_v(1,l),1)
664           massepn=ssum(iim,masse(1,l),1)
665           qpn=0.
666           do ij=1,iim
667              qpn=qpn+masse(ij,l)*q(ij,l)
668           enddo
669           qpn=(qpn+convpn)/(massepn+convmpn)
670           do ij=1,iip1
671              q(ij,l)=qpn
672           enddo
673         endif
674         
675c        convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)/apols
676c        convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)/apols
677         
678         if (pole_sud) then
679         
680           convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
681           convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
682           masseps=ssum(iim, masse(ip1jm+1,l),1)
683           qps=0.
684           do ij = ip1jm+1,ip1jmp1-1
685              qps=qps+masse(ij,l)*q(ij,l)
686           enddo
687           qps=(qps+convps)/(masseps+convmps)
688           do ij=ip1jm+1,ip1jmp1
689              q(ij,l)=qps
690           enddo
691         endif
692c.-. fin ancienne version
693
694c._. nouvelle version
695c        convpn=SSUM(iim,qbyv(1,l),1)
696c        convmpn=ssum(iim,masse_adv_v(1,l),1)
697c        oldmasse=ssum(iim,masse(1,l),1)
698c        newmasse=oldmasse+convmpn
699c        newq=(q(1,l)*oldmasse+convpn)/newmasse
700c        newmasse=newmasse/apoln
701c        DO ij = 1,iip1
702c           q(ij,l)=newq
703c           masse(ij,l)=newmasse*aire(ij)
704c        ENDDO
705c        convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
706c        convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
707c        oldmasse=ssum(iim,masse(ip1jm-iim,l),1)
708c        newmasse=oldmasse+convmps
709c        newq=(q(ip1jmp1,l)*oldmasse+convps)/newmasse
710c        newmasse=newmasse/apols
711c        DO ij = ip1jm+1,ip1jmp1
712c           q(ij,l)=newq
713c           masse(ij,l)=newmasse*aire(ij)
714c        ENDDO
715c._. fin nouvelle version
716      ENDDO
717c$OMP END DO NOWAIT
718
719      RETURN
720      END
721     
722     
723     
724      SUBROUTINE vlz_loc(q,pente_max,masse,w,ijb_x,ije_x)
725c
726c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
727c
728c    ********************************************************************
729c     Shema  d'advection " pseudo amont " .
730c    ********************************************************************
731c    q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....
732c     dq               sont des arguments de sortie pour le s-pg ....
733c
734c
735c   --------------------------------------------------------------------
736      USE Parallel
737      USE vlz_mod
738      IMPLICIT NONE
739c
740#include "dimensions.h"
741#include "paramet.h"
742#include "logic.h"
743#include "comvert.h"
744#include "comconst.h"
745c
746c
747c   Arguments:
748c   ----------
749      REAL masse(ijb_u:ije_u,llm),pente_max
750      REAL q(ijb_u:ije_u,llm)
751      REAL w(ijb_u:ije_u,llm+1)
752c
753c      Local
754c   ---------
755c
756      INTEGER i,ij,l,j,ii
757c
758      REAL newmasse
759
760      REAL dzqmax
761      REAL sigw
762
763      LOGICAL testcpu
764      SAVE testcpu
765c$OMP THREADPRIVATE(testcpu)
766      REAL temps0,temps1,temps2,temps3,temps4,temps5,second
767      SAVE temps0,temps1,temps2,temps3,temps4,temps5
768c$OMP THREADPRIVATE(temps0,temps1,temps2,temps3,temps4,temps5)
769
770      REAL      SSUM
771      EXTERNAL  SSUM
772
773      DATA testcpu/.false./
774      DATA temps0,temps1,temps2,temps3,temps4,temps5/0.,0.,0.,0.,0.,0./
775      INTEGER ijb,ije,ijb_x,ije_x
776      LOGICAL,SAVE :: first=.TRUE.
777!$OMP THREADPRIVATE(first)
778     
779
780      IF (first) THEN
781       first=.FALSE.
782      ENDIF             
783c    On oriente tout dans le sens de la pression c'est a dire dans le
784c    sens de W
785
786#ifdef BIDON
787      IF(testcpu) THEN
788         temps0=second(0.)
789      ENDIF
790#endif
791
792      ijb=ijb_x
793      ije=ije_x
794
795c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
796      DO l=2,llm
797         DO ij=ijb,ije
798            dzqw(ij,l)=q(ij,l-1)-q(ij,l)
799            adzqw(ij,l)=abs(dzqw(ij,l))
800         ENDDO
801      ENDDO
802c$OMP END DO
803
804c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
805      DO l=2,llm-1
806         DO ij=ijb,ije
807#ifdef CRAY
808            dzq(ij,l)=0.5*
809     ,      cvmgp(dzqw(ij,l)+dzqw(ij,l+1),0.,dzqw(ij,l)*dzqw(ij,l+1))
810#else
811            IF(dzqw(ij,l)*dzqw(ij,l+1).gt.0.) THEN
812                dzq(ij,l)=0.5*(dzqw(ij,l)+dzqw(ij,l+1))
813            ELSE
814                dzq(ij,l)=0.
815            ENDIF
816#endif
817            dzqmax=pente_max*min(adzqw(ij,l),adzqw(ij,l+1))
818            dzq(ij,l)=sign(min(abs(dzq(ij,l)),dzqmax),dzq(ij,l))
819         ENDDO
820      ENDDO
821c$OMP END DO NOWAIT
822
823c$OMP MASTER
824      DO ij=ijb,ije
825         dzq(ij,1)=0.
826         dzq(ij,llm)=0.
827      ENDDO
828c$OMP END MASTER
829c$OMP BARRIER
830#ifdef BIDON
831      IF(testcpu) THEN
832         temps1=temps1+second(0.)-temps0
833      ENDIF
834#endif
835c ---------------------------------------------------------------
836c   .... calcul des termes d'advection verticale  .......
837c ---------------------------------------------------------------
838
839c calcul de  - d( q   * w )/ d(sigma)    qu'on ajoute a  dq pour calculer dq
840
841c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
842       DO l = 1,llm-1
843         do  ij = ijb,ije
844          IF(w(ij,l+1).gt.0.) THEN
845             sigw=w(ij,l+1)/masse(ij,l+1)
846             wq(ij,l+1)=w(ij,l+1)*(q(ij,l+1)+0.5*(1.-sigw)*dzq(ij,l+1))
847          ELSE
848             sigw=w(ij,l+1)/masse(ij,l)
849             wq(ij,l+1)=w(ij,l+1)*(q(ij,l)-0.5*(1.+sigw)*dzq(ij,l))
850          ENDIF
851         ENDDO
852       ENDDO
853c$OMP END DO NOWAIT
854
855c$OMP MASTER
856       DO ij=ijb,ije
857          wq(ij,llm+1)=0.
858          wq(ij,1)=0.
859       ENDDO
860c$OMP END MASTER
861c$OMP BARRIER
862
863c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
864      DO l=1,llm
865         DO ij=ijb,ije
866            newmasse=masse(ij,l)+w(ij,l+1)-w(ij,l)
867            q(ij,l)=(q(ij,l)*masse(ij,l)+wq(ij,l+1)-wq(ij,l))
868     &         /newmasse
869            masse(ij,l)=newmasse
870         ENDDO
871      ENDDO
872c$OMP END DO NOWAIT
873
874
875      RETURN
876      END
877c      SUBROUTINE minmaxq(zq,qmin,qmax,comment)
878c
879c#include "dimensions.h"
880c#include "paramet.h"
881
882c      CHARACTER*(*) comment
883c      real qmin,qmax
884c      real zq(ip1jmp1,llm)
885
886c      INTEGER jadrs(ip1jmp1), jbad, k, i
887
888
889c      DO k = 1, llm
890c         jbad = 0
891c         DO i = 1, ip1jmp1
892c         IF (zq(i,k).GT.qmax .OR. zq(i,k).LT.qmin) THEN
893c            jbad = jbad + 1
894c            jadrs(jbad) = i
895c         ENDIF
896c         ENDDO
897c         IF (jbad.GT.0) THEN
898c         PRINT*, comment
899c         DO i = 1, jbad
900cc            PRINT*, "i,k,zq=", jadrs(i),k,zq(jadrs(i),k)
901c         ENDDO
902c         ENDIF
903c      ENDDO
904
905c      return
906c      end
907
908
909
910
Note: See TracBrowser for help on using the repository browser.