source: LMDZ6/branches/LMDZ_ECRad/libf/dyn3dmem/vlspltqs_loc.F @ 5428

Last change on this file since 5428 was 4482, checked in by lguez, 21 months ago

Sync latest trunk changes to branch LMDZ_ECRad

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