source: LMDZ6/trunk/libf/dyn3dmem/vlspltqs_loc.F @ 3601

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