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
Line 
1!
2!     $Id: vlspltqs_loc.F 4482 2023-03-29 13:14:27Z fhourdin $
3!
4      SUBROUTINE vlxqs_loc(q,pente_max,masse,u_m,qsat,ijb_x,ije_x,iq)
5c
6c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
7c
8c    ********************************************************************
9c     Shema  d''advection " pseudo amont " .
10c    ********************************************************************
11c
12c   --------------------------------------------------------------------
13      USE parallel_lmdz
14      USE infotrac, ONLY : nqtot,tracers, ! CRisi                 &
15     &                     min_qParent,min_qMass,min_ratio ! MVals et CRisi
16      IMPLICIT NONE
17c
18      include "dimensions.h"
19      include "paramet.h"
20c
21c
22c   Arguments:
23c   ----------
24      REAL masse(ijb_u:ije_u,llm,nqtot),pente_max
25      REAL u_m( ijb_u:ije_u,llm )
26      REAL q(ijb_u:ije_u,llm,nqtot)
27      REAL qsat(ijb_u:ije_u,llm)
28      INTEGER iq ! CRisi
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)
41      REAL Ratio(ijb_u:ije_u,llm,nqtot) ! CRisi
42      INTEGER ifils,iq2 ! CRisi
43
44
45      REAL      SSUM
46
47
48      INTEGER ijb,ije,ijb_x,ije_x
49     
50      !write(*,*) 'vlspltqs 58: entree vlxqs_loc, iq,ijb_x=',
51!     &   iq,ijb_x
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
75               dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq)
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
128               dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq)
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
175          iadvplus(ip1jm+1:ip1jmp1,l)=0
176        ENDDO
177c$OMP END DO NOWAIT
178      endif
179             
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
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),
189     ,                     u_m(ij,l))
190          zdum(ij,l)=0.5*zdum(ij,l)
191          u_mq(ij,l)=cvmgp(
192     ,                q(ij,l,iq)+zdum(ij,l)*dxq(ij,l),
193     ,                q(ij+1,l,iq)-zdum(ij,l)*dxq(ij+1,l),
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.
203c   le rapport de melange de l''air advecte est min(q_vanleer, Qsat_downwind)
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
208             zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l,iq)
209             u_mq(ij,l)=u_m(ij,l)*
210     $         min(q(ij,l,iq)+0.5*zdum(ij,l)*dxq(ij,l),qsat(ij+1,l))
211          ELSE
212             zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l,iq)
213             u_mq(ij,l)=u_m(ij,l)*
214     $         min(q(ij+1,l,iq)-0.5*zdum(ij,l)*dxq(ij+1,l),qsat(ij,l))
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
264cym ATTENTION ICI en OpenMP reduction pas forcement necessaire
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
281               !PRINT*,'vlxqs 280: niju,nl',niju,nl(l)
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
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)
297                        i=mod(i-2+iim,iim)+1
298                        ijq=(j-1)*iip1+i
299                     ENDDO
300c   ajout de la maille non completement advectee
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))
303                  ELSE
304                     ijq=ij+1
305                     i=ijq-(j-1)*iip1
306c   accumulation pour les mailles completements advectees
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)
311                        i=mod(i,iim)+1
312                        ijq=(j-1)*iip1+i
313                     ENDDO
314c   ajout de la maille non completement advectee
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))
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
335! CRisi: appel recursif de l'advection sur les fils.
336! Il faut faire ca avant d'avoir mis a jour q et masse
337      !write(*,*) 'vlspltqs 336: iq,ijb_x,nqChildren(iq)=',
338!     &     iq,ijb_x,tracers(iq)%nqChildren
339
340      do ifils=1,tracers(iq)%nqDescen
341        iq2=tracers(iq)%iqDescen(ifils)
342c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
343        DO l=1,llm
344          DO ij=ijb,ije
345            !MVals: veiller a ce qu'on n'ait pas de denominateur nul
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
348              Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
349            else
350              Ratio(ij,l,iq2)=min_ratio
351            endif
352          enddo   
353        enddo
354c$OMP END DO NOWAIT
355      enddo
356      do ifils=1,tracers(iq)%nqChildren
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
361! end CRisi
362
363      !write(*,*) 'vlspltqs 360: iq,ijb_x=',iq,ijb_x   
364
365c   calcul des tendances
366c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
367      DO l=1,llm
368         DO ij=ijb+1,ije
369            !MVals: veiller a ce qu'on n'ait pas de denominateur nul
370            new_m=max(masse(ij,l,iq)+u_m(ij-1,l)-u_m(ij,l),min_qMass)
371            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+
372     &      u_mq(ij-1,l)-u_mq(ij,l))
373     &      /new_m
374            masse(ij,l,iq)=new_m
375         ENDDO
376c   Modif Fred 22 03 96 correction d''un bug (les scopy ci-dessous)
377         DO ij=ijb+iip1-1,ije,iip1
378            q(ij-iim,l,iq)=q(ij,l,iq)
379            masse(ij-iim,l,iq)=masse(ij,l,iq)
380         ENDDO
381      ENDDO
382c$OMP END DO NOWAIT
383
384      !write(*,*) 'vlspltqs 380: iq,ijb_x=',iq,ijb_x
385
386! retablir les fils en rapport de melange par rapport a l'air:
387      do ifils=1,tracers(iq)%nqDescen
388        iq2=tracers(iq)%iqDescen(ifils)
389c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
390        DO l=1,llm
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
395            q(ij-iim,l,iq2)=q(ij,l,iq2)
396          enddo ! DO ij=ijb+iip1-1,ije,iip1
397        enddo
398c$OMP END DO NOWAIT
399      enddo
400
401      !write(*,*) 'vlspltqs 399: iq,ijb_x=',iq,ijb_x
402
403c     CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
404c     CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1,iq),iip1,masse(iip2,1,iq),iip1)
405
406
407      RETURN
408      END
409      SUBROUTINE vlyqs_loc(q,pente_max,masse,masse_adv_v,qsat,iq)
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 ....
417c     qsat                est   un argument de sortie pour le s-pg ....
418c
419c
420c   --------------------------------------------------------------------
421      USE parallel_lmdz
422      USE infotrac, ONLY : nqtot,tracers, ! CRisi                 &
423     &                     min_qParent,min_qMass,min_ratio ! MVals et CRisi
424      USE comconst_mod, ONLY: pi
425      IMPLICIT NONE
426c
427      include "dimensions.h"
428      include "paramet.h"
429      include "comgeom.h"
430      include "iniprint.h" 
431c
432c
433c   Arguments:
434c   ----------
435      REAL masse(ijb_u:ije_u,llm,nqtot),pente_max
436      REAL masse_adv_v( ijb_v:ije_v,llm)
437      REAL q(ijb_u:ije_u,llm,nqtot)
438      REAL qsat(ijb_u:ije_u,llm)
439      INTEGER iq ! CRisi
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)
449      REAL qbyv(ijb_v:ije_v,llm,nqtot)
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
465      REAL Ratio(ijb_u:ije_u,llm,nqtot) ! CRisi
466      INTEGER ifils,iq2 ! CRisi
467
468      REAL      SSUM
469
470      DATA first/.true./
471      INTEGER ijb,ije
472      INTEGER ijbm,ijem
473
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
481        !write(*,*) 'vlyqs 480: ij,l,iq,ijb,q(ij,l,:)=',
482!     &             ij,l,iq,ijb,q(ij,l,:)
483      endif 
484
485      IF(first) THEN
486         PRINT*,'Shema  Amont nouveau  appele dans  Vanleer   '
487         PRINT*,'vlyqs_loc, iq=',iq
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
518          airescb(i) = aire(i+ iip1) * q(i+ iip1,l,iq)
519        ENDDO
520        qpns   = SSUM( iim,  airescb ,1 ) / airej2
521      endif
522     
523      if (pole_sud) then
524        DO i = 1, iim
525          airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l,iq)
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
539         dyqv(ij)=q(ij,l,iq)-q(ij+iip1,l,iq)
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
561           dyq(ij,l)=qpns-q(ij+iip1,l,iq)
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
586         
587      ENDIF
588     
589      IF (pole_sud) THEN
590
591        DO ij=1,iip1
592           dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l,iq)-qpsn
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       
607c   calcul des pentes limites aux poles       
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
618       
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)
628C     appn=abs(dyq(1)/dyqv(iip1+1))
629C     PRINT*,dyq(ip1jm+1)
630C     PRINT*,dyqv(ip1jm-iip1+1)
631C     apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
632C     DO ij=2,iim
633C        appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
634C        apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
635C     ENDDO
636C     appn=min(pente_max/appn,1.)
637C     apps=min(pente_max/apps,1.)
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.)
643C    &   appn=0.
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.)
646C    &   apps=0.
647C
648C   limitation des pentes aux poles
649C     DO ij=1,iip1
650C        dyq(ij)=appn*dyq(ij)
651C        dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
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
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)))
718         ELSE
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)) )
721         ENDIF
722          qbyv(ij,l,iq) = masse_adv_v(ij,l)*qbyv(ij,l,iq)
723       ENDDO
724      ENDDO
725c$OMP END DO NOWAIT
726
727! CRisi: appel recursif de l'advection sur les fils.
728! Il faut faire ca avant d'avoir mis a jour q et masse
729!     write(*,*)'vlyqs 689: iq,nqChildren(iq)=',iq,
730!    &             tracers(iq)%nqChildren
731     
732      ijb=ij_begin-2*iip1
733      ije=ij_end+2*iip1
734      ijbm=ij_begin-iip1
735      ijem=ij_end+iip1
736      if (pole_nord) ijb=ij_begin
737      if (pole_sud)  ije=ij_end 
738      if (pole_nord) ijbm=ij_begin
739      if (pole_sud)  ijem=ij_end
740
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
744      do ifils=1,tracers(iq)%nqDescen
745        iq2=tracers(iq)%iqDescen(ifils)
746c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
747        DO l=1,llm
748          ! modif des bornes: CRisi 16 nov 2020
749          ! d'abord masse avec bornes corrigees
750          DO ij=ijbm,ijem
751            !MVals: veiller a ce qu'on n'ait pas de denominateur nul
752            masse(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),min_qMass)
753          enddo !DO ij=ijbm,ijem
754
755          ! ensuite Ratio avec anciennes bornes
756          DO ij=ijb,ije
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)
759            if (q(ij,l,iq).gt.min_qParent) then ! modif 13 nov 2020
760              Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
761            else
762              Ratio(ij,l,iq2)=min_ratio   
763            endif
764          enddo !DO ij=ijbm,ijem
765        enddo !DO l=1,llm
766c$OMP END DO NOWAIT
767      enddo
768      do ifils=1,tracers(iq)%nqChildren
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
773
774       
775! end CRisi
776
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
785            newmasse=masse(ij,l,iq)
786     &      +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
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
790         ENDDO
791c.-. ancienne version
792
793         IF (pole_nord) THEN
794
795           convpn=SSUM(iim,qbyv(1,l,iq),1)/apoln
796           convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln
797           DO ij = 1,iip1
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))/
800     &                 newmasse
801              masse(ij,l,iq)=newmasse
802           ENDDO
803         
804         ENDIF
805         
806         IF (pole_sud) THEN
807         
808           convps  = -SSUM(iim,qbyv(ip1jm-iim,l,iq),iq,1)/apols
809           convmps = -SSUM(iim,masse_adv_v(ip1jm-iim,l),1)/apols
810           DO ij = ip1jm+1,ip1jmp1
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))/
813     &                 newmasse
814              masse(ij,l,iq)=newmasse
815           ENDDO
816         
817         ENDIF
818c.-. fin ancienne version
819
820c._. nouvelle version
821c        convpn=SSUM(iim,qbyv(1,l,iq),1)
822c        convmpn=ssum(iim,masse_adv_v(1,l),1)
823c        oldmasse=ssum(iim,masse(1,l,iq),1)
824c        newmasse=oldmasse+convmpn
825c        newq=(q(1,l,iq)*oldmasse+convpn)/newmasse
826c        newmasse=newmasse/apoln
827c        DO ij = 1,iip1
828c           q(ij,l,iq)=newq
829c           masse(ij,l,iq)=newmasse*aire(ij)
830c        ENDDO
831c        convps=-SSUM(iim,qbyv(ip1jm-iim,l,iq),1)
832c        convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
833c        oldmasse=ssum(iim,masse(ip1jm-iim,l,iq),1)
834c        newmasse=oldmasse+convmps
835c        newq=(q(ip1jmp1,l,iq)*oldmasse+convps)/newmasse
836c        newmasse=newmasse/apols
837c        DO ij = ip1jm+1,ip1jmp1
838c           q(ij,l,iq)=newq
839c           masse(ij,l,iq)=newmasse*aire(ij)
840c        ENDDO
841c._. fin nouvelle version
842      ENDDO
843c$OMP END DO NOWAIT
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 
851      do ifils=1,tracers(iq)%nqDescen
852        iq2=tracers(iq)%iqDescen(ifils)
853c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
854        DO l=1,llm
855          DO ij=ijb,ije
856            q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)           
857          enddo
858        enddo
859c$OMP END DO NOWAIT
860      enddo
861
862
863      RETURN
864      END
Note: See TracBrowser for help on using the repository browser.