source: LMDZ5/trunk/libf/dyn3dmem/vlsplt_loc.F @ 2740

Last change on this file since 2740 was 2603, checked in by Ehouarn Millour, 8 years ago

Cleanup in the dynamics: turn logic.h into module logic_mod.F90
EM

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