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

Last change on this file since 3836 was 3800, checked in by Laurent Fairhead, 4 years ago

Modifications nécessaires pour les isotopes
CRisi

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