source: LMDZ6/branches/Ocean_skin/libf/dyn3d/vlspltqs.F @ 4738

Last change on this file since 4738 was 4368, checked in by lguez, 2 years ago

Sync latest trunk changes to Ocean_skin

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