source: LMDZ5/trunk/libf/dyn3dmem/vlspltqs_loc.F @ 2587

Last change on this file since 2587 was 2286, checked in by Ehouarn Millour, 9 years ago
  • fix in dyn3d the array out of bound issue (in qminimum) that was only corrected in dyn3dmem in rev 2285.
  • comment out many invasive debug writes.

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