source: LMDZ4/branches/pre_V3/libf/phylmd/thermcell.F @ 789

Last change on this file since 789 was 542, checked in by lmdzadmin, 20 years ago

Oubli
LF

  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 40.8 KB
Line 
1      SUBROUTINE thermcell(ngrid,nlay,ptimestep
2     s                  ,pplay,pplev,pphi
3     s                  ,pu,pv,pt,po
4     s                  ,pduadj,pdvadj,pdtadj,pdoadj
5     s                  ,fm0,entr0
6c    s                  ,pu_therm,pv_therm
7     s                  ,r_aspect,l_mix,w2di,tho)
8
9      IMPLICIT NONE
10
11c=======================================================================
12c
13c   Calcul du transport verticale dans la couche limite en presence
14c   de "thermiques" explicitement representes
15c
16c   Réécriture à partir d'un listing papier à Habas, le 14/02/00
17c
18c   le thermique est supposé homogène et dissipé par mélange avec
19c   son environnement. la longueur l_mix contrôle l'efficacité du
20c   mélange
21c
22c   Le calcul du transport des différentes espèces se fait en prenant
23c   en compte:
24c     1. un flux de masse montant
25c     2. un flux de masse descendant
26c     3. un entrainement
27c     4. un detrainement
28c
29c=======================================================================
30
31c-----------------------------------------------------------------------
32c   declarations:
33c   -------------
34
35#include "dimensions.h"
36#include "dimphy.h"
37#include "YOMCST.h"
38
39c   arguments:
40c   ----------
41
42      INTEGER ngrid,nlay,w2di,tho
43      real ptimestep,l_mix,r_aspect
44      REAL pt(ngrid,nlay),pdtadj(ngrid,nlay)
45      REAL pu(ngrid,nlay),pduadj(ngrid,nlay)
46      REAL pv(ngrid,nlay),pdvadj(ngrid,nlay)
47      REAL po(ngrid,nlay),pdoadj(ngrid,nlay)
48      REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
49      real pphi(ngrid,nlay)
50
51      integer idetr
52      save idetr
53      data idetr/3/
54
55c   local:
56c   ------
57
58      INTEGER ig,k,l,lmaxa(klon),lmix(klon)
59      real zsortie1d(klon)
60c CR: on remplace lmax(klon,klev+1)
61      INTEGER lmax(klon),lmin(klon),lentr(klon)
62      real linter(klon)
63      real zmix(klon), fracazmix(klon)
64c RC
65      real zmax(klon),zw,zz,zw2(klon,klev+1),ztva(klon,klev),zzz
66
67      real zlev(klon,klev+1),zlay(klon,klev)
68      REAL zh(klon,klev),zdhadj(klon,klev)
69      REAL ztv(klon,klev)
70      real zu(klon,klev),zv(klon,klev),zo(klon,klev)
71      REAL wh(klon,klev+1)
72      real wu(klon,klev+1),wv(klon,klev+1),wo(klon,klev+1)
73      real zla(klon,klev+1)
74      real zwa(klon,klev+1)
75      real zld(klon,klev+1)
76      real zwd(klon,klev+1)
77      real zsortie(klon,klev)
78      real zva(klon,klev)
79      real zua(klon,klev)
80      real zoa(klon,klev)
81
82      real zha(klon,klev)
83      real wa_moy(klon,klev+1)
84      real fraca(klon,klev+1)
85      real fracc(klon,klev+1)
86      real zf,zf2
87      real thetath2(klon,klev),wth2(klon,klev)
88      common/comtherm/thetath2,wth2
89
90      real count_time
91      integer isplit,nsplit,ialt
92      parameter (nsplit=10)
93      data isplit/0/
94      save isplit
95
96      logical sorties
97      real rho(klon,klev),rhobarz(klon,klev+1),masse(klon,klev)
98      real zpspsk(klon,klev)
99
100c     real wmax(klon,klev),wmaxa(klon)
101      real wmax(klon),wmaxa(klon)
102      real wa(klon,klev,klev+1)
103      real wd(klon,klev+1)
104      real larg_part(klon,klev,klev+1)
105      real fracd(klon,klev+1)
106      real xxx(klon,klev+1)
107      real larg_cons(klon,klev+1)
108      real larg_detr(klon,klev+1)
109      real fm0(klon,klev+1),entr0(klon,klev),detr(klon,klev)
110      real pu_therm(klon,klev),pv_therm(klon,klev)
111      real fm(klon,klev+1),entr(klon,klev)
112      real fmc(klon,klev+1)
113
114cCR:nouvelles variables
115      real f_star(klon,klev+1),entr_star(klon,klev)
116      real entr_star_tot(klon),entr_star2(klon)
117      real f(klon), f0(klon)
118      real zlevinter(klon)
119      logical first
120      data first /.false./
121      save first
122cRC
123
124      character*2 str2
125      character*10 str10
126
127      LOGICAL vtest(klon),down
128
129      EXTERNAL SCOPY
130
131      integer ncorrec,ll
132      save ncorrec
133      data ncorrec/0/
134     
135c
136c-----------------------------------------------------------------------
137c   initialisation:
138c   ---------------
139c
140       sorties=.true.
141      IF(ngrid.NE.klon) THEN
142         PRINT*
143         PRINT*,'STOP dans convadj'
144         PRINT*,'ngrid    =',ngrid
145         PRINT*,'klon  =',klon
146      ENDIF
147c
148c-----------------------------------------------------------------------
149c   incrementation eventuelle de tendances precedentes:
150c   ---------------------------------------------------
151
152       print*,'0 OK convect8'
153
154      DO 1010 l=1,nlay
155         DO 1015 ig=1,ngrid
156            zpspsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**RKAPPA
157            zh(ig,l)=pt(ig,l)/zpspsk(ig,l)
158            zu(ig,l)=pu(ig,l)
159            zv(ig,l)=pv(ig,l)
160            zo(ig,l)=po(ig,l)
161            ztv(ig,l)=zh(ig,l)*(1.+0.61*zo(ig,l))
1621015     CONTINUE
1631010  CONTINUE
164
165       print*,'1 OK convect8'
166c                       --------------------
167c
168c
169c                       + + + + + + + + + + +
170c
171c
172c  wa, fraca, wd, fracd --------------------   zlev(2), rhobarz
173c  wh,wt,wo ...
174c
175c                       + + + + + + + + + + +  zh,zu,zv,zo,rho
176c
177c
178c                       --------------------   zlev(1)
179c                       \\\\\\\\\\\\\\\\\\\\
180c
181c
182
183c-----------------------------------------------------------------------
184c   Calcul des altitudes des couches
185c-----------------------------------------------------------------------
186
187      do l=2,nlay
188         do ig=1,ngrid
189            zlev(ig,l)=0.5*(pphi(ig,l)+pphi(ig,l-1))/RG
190         enddo
191      enddo
192      do ig=1,ngrid
193         zlev(ig,1)=0.
194         zlev(ig,nlay+1)=(2.*pphi(ig,klev)-pphi(ig,klev-1))/RG
195      enddo
196      do l=1,nlay
197         do ig=1,ngrid
198            zlay(ig,l)=pphi(ig,l)/RG
199         enddo
200      enddo
201
202c      print*,'2 OK convect8'
203c-----------------------------------------------------------------------
204c   Calcul des densites
205c-----------------------------------------------------------------------
206
207      do l=1,nlay
208         do ig=1,ngrid
209            rho(ig,l)=pplay(ig,l)/(zpspsk(ig,l)*RD*zh(ig,l))
210         enddo
211      enddo
212
213      do l=2,nlay
214         do ig=1,ngrid
215            rhobarz(ig,l)=0.5*(rho(ig,l)+rho(ig,l-1))
216         enddo
217      enddo
218
219      do k=1,nlay
220         do l=1,nlay+1
221            do ig=1,ngrid
222               wa(ig,k,l)=0.
223            enddo
224         enddo
225      enddo
226
227c      print*,'3 OK convect8'
228c------------------------------------------------------------------
229c   Calcul de w2, quarre de w a partir de la cape
230c   a partir de w2, on calcule wa, vitesse de l'ascendance
231c
232c   ATTENTION: Dans cette version, pour cause d'economie de memoire,
233c   w2 est stoke dans wa
234c
235c   ATTENTION: dans convect8, on n'utilise le calcule des wa
236c   independants par couches que pour calculer l'entrainement
237c   a la base et la hauteur max de l'ascendance.
238c
239c   Indicages:
240c   l'ascendance provenant du niveau k traverse l'interface l avec
241c   une vitesse wa(k,l).
242c
243c                       --------------------
244c
245c                       + + + + + + + + + +
246c
247c  wa(k,l)   ----       --------------------    l
248c             /\
249c            /||\       + + + + + + + + + +
250c             ||
251c             ||        --------------------
252c             ||
253c             ||        + + + + + + + + + +
254c             ||
255c             ||        --------------------
256c             ||__
257c             |___      + + + + + + + + + +     k
258c
259c                       --------------------
260c
261c
262c
263c------------------------------------------------------------------
264
265cCR: ponderation entrainement des couches instables
266cdef des entr_star tels que entr=f*entr_star     
267      do l=1,klev
268         do ig=1,ngrid
269            entr_star(ig,l)=0.
270         enddo
271      enddo
272c determination de la longueur de la couche d entrainement
273      do ig=1,ngrid
274         lentr(ig)=1
275      enddo
276
277con ne considere que les premieres couches instables
278      do k=nlay-2,1,-1
279         do ig=1,ngrid
280            if (ztv(ig,k).gt.ztv(ig,k+1).and.
281     s          ztv(ig,k+1).le.ztv(ig,k+2)) then
282               lentr(ig)=k
283            endif
284          enddo
285      enddo
286   
287c determination du lmin: couche d ou provient le thermique
288      do ig=1,ngrid
289         lmin(ig)=1
290      enddo
291      do ig=1,ngrid
292         do l=nlay,2,-1
293            if (ztv(ig,l-1).gt.ztv(ig,l)) then
294               lmin(ig)=l-1
295            endif
296         enddo
297      enddo
298c
299c definition de l'entrainement des couches
300      do l=1,klev-1
301         do ig=1,ngrid
302            if (ztv(ig,l).gt.ztv(ig,l+1).and.
303     s          l.ge.lmin(ig).and.l.le.lentr(ig)) then
304                 entr_star(ig,l)=(ztv(ig,l)-ztv(ig,l+1))*
305     s                           (zlev(ig,l+1)-zlev(ig,l))
306            endif
307         enddo
308      enddo
309c pas de thermique si couches 1->5 stables
310      do ig=1,ngrid
311         if (lmin(ig).gt.5) then
312            do l=1,klev
313               entr_star(ig,l)=0.
314            enddo
315         endif
316      enddo
317c calcul de l entrainement total
318      do ig=1,ngrid
319         entr_star_tot(ig)=0.
320      enddo
321      do ig=1,ngrid
322         do k=1,klev
323            entr_star_tot(ig)=entr_star_tot(ig)+entr_star(ig,k)
324         enddo
325      enddo
326c
327      print*,'fin calcul entr_star'
328      do k=1,klev
329         do ig=1,ngrid
330            ztva(ig,k)=ztv(ig,k)
331         enddo
332      enddo
333cRC
334c      print*,'7 OK convect8'
335      do k=1,klev+1
336         do ig=1,ngrid
337            zw2(ig,k)=0.
338            fmc(ig,k)=0.
339cCR
340            f_star(ig,k)=0.
341cRC
342            larg_cons(ig,k)=0.
343            larg_detr(ig,k)=0.
344            wa_moy(ig,k)=0.
345         enddo
346      enddo
347
348c      print*,'8 OK convect8'
349      do ig=1,ngrid
350         linter(ig)=1.
351         lmaxa(ig)=1
352         lmix(ig)=1
353         wmaxa(ig)=0.
354      enddo
355
356cCR:
357      do l=1,nlay-2
358         do ig=1,ngrid
359            if (ztv(ig,l).gt.ztv(ig,l+1)
360     s         .and.entr_star(ig,l).gt.1.e-10
361     s         .and.zw2(ig,l).lt.1e-10) then
362               f_star(ig,l+1)=entr_star(ig,l)
363ctest:calcul de dteta
364               zw2(ig,l+1)=2.*RG*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1)
365     s                     *(zlev(ig,l+1)-zlev(ig,l))
366     s                     *0.4*pphi(ig,l)/(pphi(ig,l+1)-pphi(ig,l))
367               larg_detr(ig,l)=0.
368            else if ((zw2(ig,l).ge.1e-10).and.
369     s               (f_star(ig,l)+entr_star(ig,l).gt.1.e-10)) then
370               f_star(ig,l+1)=f_star(ig,l)+entr_star(ig,l)
371               ztva(ig,l)=(f_star(ig,l)*ztva(ig,l-1)+entr_star(ig,l)
372     s                    *ztv(ig,l))/f_star(ig,l+1)
373               zw2(ig,l+1)=zw2(ig,l)*(f_star(ig,l)/f_star(ig,l+1))**2+
374     s                     2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
375     s                     *(zlev(ig,l+1)-zlev(ig,l))
376            endif
377c determination de zmax continu par interpolation lineaire
378            if (zw2(ig,l+1).lt.0.) then
379ctest
380               if (abs(zw2(ig,l+1)-zw2(ig,l)).lt.1e-10) then
381                  print*,'pb linter'
382               endif
383               linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))
384     s           -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
385               zw2(ig,l+1)=0.
386               lmaxa(ig)=l
387            else
388               if (zw2(ig,l+1).lt.0.) then
389                  print*,'pb1 zw2<0'
390               endif
391               wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
392            endif
393            if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
394c   lmix est le niveau de la couche ou w (wa_moy) est maximum
395               lmix(ig)=l+1
396               wmaxa(ig)=wa_moy(ig,l+1)
397            endif
398         enddo
399      enddo
400      print*,'fin calcul zw2'
401c
402c Calcul de la couche correspondant a la hauteur du thermique
403      do ig=1,ngrid
404         lmax(ig)=lentr(ig)
405      enddo
406      do ig=1,ngrid
407         do l=nlay,lentr(ig)+1,-1
408            if (zw2(ig,l).le.1.e-10) then
409               lmax(ig)=l-1
410            endif
411         enddo
412      enddo
413c pas de thermique si couches 1->5 stables
414      do ig=1,ngrid
415         if (lmin(ig).gt.5) then
416            lmax(ig)=1
417            lmin(ig)=1
418         endif
419      enddo
420c   
421c Determination de zw2 max
422      do ig=1,ngrid
423         wmax(ig)=0.
424      enddo
425
426      do l=1,nlay
427         do ig=1,ngrid
428            if (l.le.lmax(ig)) then
429                if (zw2(ig,l).lt.0.)then
430                  print*,'pb2 zw2<0'
431                endif
432                zw2(ig,l)=sqrt(zw2(ig,l))
433                wmax(ig)=max(wmax(ig),zw2(ig,l))
434            else
435                 zw2(ig,l)=0.
436            endif
437          enddo
438      enddo
439
440c   Longueur caracteristique correspondant a la hauteur des thermiques.
441      do  ig=1,ngrid
442         zmax(ig)=0.
443         zlevinter(ig)=zlev(ig,1)
444      enddo
445      do  ig=1,ngrid
446c calcul de zlevinter
447          zlevinter(ig)=(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*
448     s    linter(ig)+zlev(ig,lmax(ig))-lmax(ig)*(zlev(ig,lmax(ig)+1)
449     s    -zlev(ig,lmax(ig)))
450       zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,lmin(ig)))
451      enddo
452
453      print*,'avant fermeture'
454c Fermeture,determination de f
455      do ig=1,ngrid
456         entr_star2(ig)=0.
457      enddo
458      do ig=1,ngrid
459         if (entr_star_tot(ig).LT.1.e-10) then
460            f(ig)=0.
461         else
462             do k=lmin(ig),lentr(ig)
463                entr_star2(ig)=entr_star2(ig)+entr_star(ig,k)**2
464     s                    /(rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k)))
465             enddo
466c Nouvelle fermeture
467             f(ig)=wmax(ig)/(max(500.,zmax(ig))*r_aspect
468     s             *entr_star2(ig))*entr_star_tot(ig)
469ctest
470c             if (first) then
471c             f(ig)=f(ig)+(f0(ig)-f(ig))*exp(-ptimestep/zmax(ig)
472c     s             *wmax(ig))
473c             endif
474         endif
475c         f0(ig)=f(ig)
476c         first=.true.
477      enddo
478      print*,'apres fermeture'
479
480c Calcul de l'entrainement
481       do k=1,klev
482         do ig=1,ngrid
483            entr(ig,k)=f(ig)*entr_star(ig,k)
484         enddo
485      enddo
486c Calcul des flux
487      do ig=1,ngrid
488         do l=1,lmax(ig)-1
489            fmc(ig,l+1)=fmc(ig,l)+entr(ig,l)
490         enddo
491      enddo
492
493cRC
494
495
496c      print*,'9 OK convect8'
497c     print*,'WA1 ',wa_moy
498
499c   determination de l'indice du debut de la mixed layer ou w decroit
500
501c   calcul de la largeur de chaque ascendance dans le cas conservatif.
502c   dans ce cas simple, on suppose que la largeur de l'ascendance provenant
503c   d'une couche est égale à la hauteur de la couche alimentante.
504c   La vitesse maximale dans l'ascendance est aussi prise comme estimation
505c   de la vitesse d'entrainement horizontal dans la couche alimentante.
506
507      do l=2,nlay
508         do ig=1,ngrid
509            if (l.le.lmaxa(ig)) then
510               zw=max(wa_moy(ig,l),1.e-10)
511               larg_cons(ig,l)=zmax(ig)*r_aspect
512     s         *fmc(ig,l)/(rhobarz(ig,l)*zw)
513            endif
514         enddo
515      enddo
516
517      do l=2,nlay
518         do ig=1,ngrid
519            if (l.le.lmaxa(ig)) then
520c              if (idetr.eq.0) then
521c  cette option est finalement en dur.
522                  if ((l_mix*zlev(ig,l)).lt.0.)then
523                   print*,'pb l_mix*zlev<0'
524                  endif
525                  larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
526c              else if (idetr.eq.1) then
527c                 larg_detr(ig,l)=larg_cons(ig,l)
528c    s            *sqrt(l_mix*zlev(ig,l))/larg_cons(ig,lmix(ig))
529c              else if (idetr.eq.2) then
530c                 larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
531c    s            *sqrt(wa_moy(ig,l))
532c              else if (idetr.eq.4) then
533c                 larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
534c    s            *wa_moy(ig,l)
535c              endif
536            endif
537         enddo
538       enddo
539
540c      print*,'10 OK convect8'
541c     print*,'WA2 ',wa_moy
542c   calcul de la fraction de la maille concernée par l'ascendance en tenant
543c   compte de l'epluchage du thermique.
544c
545cCR def de  zmix continu (profil parabolique des vitesses)
546      do ig=1,ngrid
547           if (lmix(ig).gt.1.) then
548c test
549              if (((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
550     s        *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1)))
551     s        -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
552     s        *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))).gt.1e-10)
553     s        then
554c             
555            zmix(ig)=((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
556     s        *((zlev(ig,lmix(ig)))**2-(zlev(ig,lmix(ig)+1))**2)
557     s        -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
558     s        *((zlev(ig,lmix(ig)-1))**2-(zlev(ig,lmix(ig)))**2))
559     s        /(2.*((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
560     s        *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1)))
561     s        -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
562     s        *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))))
563            else
564            zmix(ig)=zlev(ig,lmix(ig))
565            print*,'pb zmix'
566            endif
567         else
568         zmix(ig)=0.
569         endif
570ctest
571         if ((zmax(ig)-zmix(ig)).lt.0.) then
572            zmix(ig)=0.99*zmax(ig)
573c            print*,'pb zmix>zmax'
574         endif
575      enddo
576c
577c calcul du nouveau lmix correspondant
578      do ig=1,ngrid
579         do l=1,klev
580            if (zmix(ig).ge.zlev(ig,l).and.
581     s          zmix(ig).lt.zlev(ig,l+1)) then
582              lmix(ig)=l
583             endif
584          enddo
585      enddo
586c
587      do l=2,nlay
588         do ig=1,ngrid
589            if(larg_cons(ig,l).gt.1.) then
590c     print*,ig,l,lmix(ig),lmaxa(ig),larg_cons(ig,l),'  KKK'
591               fraca(ig,l)=(larg_cons(ig,l)-larg_detr(ig,l))
592     s            /(r_aspect*zmax(ig))
593c test
594               fraca(ig,l)=max(fraca(ig,l),0.)
595               fraca(ig,l)=min(fraca(ig,l),0.5)
596               fracd(ig,l)=1.-fraca(ig,l)
597               fracc(ig,l)=larg_cons(ig,l)/(r_aspect*zmax(ig))
598            else
599c              wa_moy(ig,l)=0.
600               fraca(ig,l)=0.
601               fracc(ig,l)=0.
602               fracd(ig,l)=1.
603            endif
604         enddo
605      enddo                 
606cCR: calcul de fracazmix
607       do ig=1,ngrid
608          fracazmix(ig)=(fraca(ig,lmix(ig)+1)-fraca(ig,lmix(ig)))/
609     s     (zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))*zmix(ig)
610     s    +fraca(ig,lmix(ig))-zlev(ig,lmix(ig))*(fraca(ig,lmix(ig)+1)
611     s    -fraca(ig,lmix(ig)))/(zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))
612       enddo
613c
614       do l=2,nlay
615          do ig=1,ngrid
616             if(larg_cons(ig,l).gt.1.) then
617               if (l.gt.lmix(ig)) then
618ctest
619                 if (zmax(ig)-zmix(ig).lt.1.e-10) then
620c                   print*,'pb xxx'
621                   xxx(ig,l)=(lmaxa(ig)+1.-l)/(lmaxa(ig)+1.-lmix(ig))
622                 else
623                 xxx(ig,l)=(zmax(ig)-zlev(ig,l))/(zmax(ig)-zmix(ig))
624                 endif
625           if (idetr.eq.0) then
626               fraca(ig,l)=fracazmix(ig)
627           else if (idetr.eq.1) then
628               fraca(ig,l)=fracazmix(ig)*xxx(ig,l)
629           else if (idetr.eq.2) then
630               fraca(ig,l)=fracazmix(ig)*(1.-(1.-xxx(ig,l))**2)
631           else
632               fraca(ig,l)=fracazmix(ig)*xxx(ig,l)**2
633           endif
634c     print*,ig,l,lmix(ig),lmaxa(ig),xxx(ig,l),'LLLLLLL'
635               fraca(ig,l)=max(fraca(ig,l),0.)
636               fraca(ig,l)=min(fraca(ig,l),0.5)
637               fracd(ig,l)=1.-fraca(ig,l)
638               fracc(ig,l)=larg_cons(ig,l)/(r_aspect*zmax(ig))
639             endif
640            endif
641         enddo
642      enddo
643     
644      print*,'fin calcul fraca'
645c      print*,'11 OK convect8'
646c     print*,'Ea3 ',wa_moy
647c------------------------------------------------------------------
648c   Calcul de fracd, wd
649c   somme wa - wd = 0
650c------------------------------------------------------------------
651
652
653      do ig=1,ngrid
654         fm(ig,1)=0.
655         fm(ig,nlay+1)=0.
656      enddo
657
658      do l=2,nlay
659           do ig=1,ngrid
660              fm(ig,l)=fraca(ig,l)*wa_moy(ig,l)*rhobarz(ig,l)
661cCR:test
662              if (entr(ig,l-1).lt.1e-10.and.fm(ig,l).gt.fm(ig,l-1)
663     s            .and.l.gt.lmix(ig)) then
664                 fm(ig,l)=fm(ig,l-1)
665c                 write(1,*)'ajustement fm, l',l
666              endif
667c              write(1,*)'ig,l,fm(ig,l)',ig,l,fm(ig,l)
668cRC
669           enddo
670         do ig=1,ngrid
671            if(fracd(ig,l).lt.0.1) then
672               stop'fracd trop petit'
673            else
674c    vitesse descendante "diagnostique"
675               wd(ig,l)=fm(ig,l)/(fracd(ig,l)*rhobarz(ig,l))
676            endif
677         enddo
678      enddo
679
680      do l=1,nlay
681         do ig=1,ngrid
682c           masse(ig,l)=rho(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
683            masse(ig,l)=(pplev(ig,l)-pplev(ig,l+1))/RG
684         enddo
685      enddo
686
687      print*,'12 OK convect8'
688c     print*,'WA4 ',wa_moy
689cc------------------------------------------------------------------
690c   calcul du transport vertical
691c------------------------------------------------------------------
692
693      go to 4444
694c     print*,'XXXXXXXXXXXXXXX ptimestep= ',ptimestep
695      do l=2,nlay-1
696         do ig=1,ngrid
697            if(fm(ig,l+1)*ptimestep.gt.masse(ig,l)
698     s      .and.fm(ig,l+1)*ptimestep.gt.masse(ig,l+1)) then
699c     print*,'WARN!!! FM>M ig=',ig,' l=',l,'  FM='
700c    s         ,fm(ig,l+1)*ptimestep
701c    s         ,'   M=',masse(ig,l),masse(ig,l+1)
702             endif
703         enddo
704      enddo
705
706      do l=1,nlay
707         do ig=1,ngrid
708            if(entr(ig,l)*ptimestep.gt.masse(ig,l)) then
709c     print*,'WARN!!! E>M ig=',ig,' l=',l,'  E=='
710c    s         ,entr(ig,l)*ptimestep
711c    s         ,'   M=',masse(ig,l)
712             endif
713         enddo
714      enddo
715
716      do l=1,nlay
717         do ig=1,ngrid
718            if(.not.fm(ig,l).ge.0..or..not.fm(ig,l).le.10.) then
719c     print*,'WARN!!! fm exagere ig=',ig,'   l=',l
720c    s         ,'   FM=',fm(ig,l)
721            endif
722            if(.not.masse(ig,l).ge.1.e-10
723     s         .or..not.masse(ig,l).le.1.e4) then
724c     print*,'WARN!!! masse exagere ig=',ig,'   l=',l
725c    s         ,'   M=',masse(ig,l)
726c     print*,'rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)',
727c    s                 rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)
728c     print*,'zlev(ig,l+1),zlev(ig,l)'
729c    s                ,zlev(ig,l+1),zlev(ig,l)
730c     print*,'pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)'
731c    s                ,pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)
732            endif
733            if(.not.entr(ig,l).ge.0..or..not.entr(ig,l).le.10.) then
734c     print*,'WARN!!! entr exagere ig=',ig,'   l=',l
735c    s         ,'   E=',entr(ig,l)
736            endif
737         enddo
738      enddo
739
7404444   continue
741
742cCR:redefinition du entr
743       do l=1,nlay
744         do ig=1,ngrid
745            detr(ig,l)=fm(ig,l)+entr(ig,l)-fm(ig,l+1)
746            if (detr(ig,l).lt.0.) then
747                entr(ig,l)=entr(ig,l)-detr(ig,l)
748                detr(ig,l)=0.
749c     print*,'WARNING !!! detrainement negatif ',ig,l
750            endif
751         enddo
752      enddo
753cRC
754      if (w2di.eq.1) then
755         fm0=fm0+ptimestep*(fm-fm0)/float(tho)
756         entr0=entr0+ptimestep*(entr-entr0)/float(tho)
757      else
758         fm0=fm
759         entr0=entr
760      endif
761
762      if (1.eq.1) then
763         call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
764     .    ,zh,zdhadj,zha)
765         call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
766     .    ,zo,pdoadj,zoa)
767      else
768         call dqthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca
769     .    ,zh,zdhadj,zha)
770         call dqthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca
771     .    ,zo,pdoadj,zoa)
772      endif
773
774      if (1.eq.0) then
775         call dvthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse
776     .    ,fraca,zmax
777     .    ,zu,zv,pduadj,pdvadj,zua,zva)
778      else
779         call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
780     .    ,zu,pduadj,zua)
781         call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
782     .    ,zv,pdvadj,zva)
783      endif
784
785      do l=1,nlay
786         do ig=1,ngrid
787            zf=0.5*(fracc(ig,l)+fracc(ig,l+1))
788            zf2=zf/(1.-zf)
789            thetath2(ig,l)=zf2*(zha(ig,l)-zh(ig,l))**2
790            wth2(ig,l)=zf2*(0.5*(wa_moy(ig,l)+wa_moy(ig,l+1)))**2
791         enddo
792      enddo
793
794
795
796c     print*,'13 OK convect8'
797c     print*,'WA5 ',wa_moy
798      do l=1,nlay
799         do ig=1,ngrid
800            pdtadj(ig,l)=zdhadj(ig,l)*zpspsk(ig,l)
801         enddo
802      enddo
803
804
805c     do l=1,nlay
806c        do ig=1,ngrid
807c           if(abs(pdtadj(ig,l))*86400..gt.500.) then
808c     print*,'WARN!!! ig=',ig,'  l=',l
809c    s         ,'   pdtadj=',pdtadj(ig,l)
810c           endif
811c           if(abs(pdoadj(ig,l))*86400..gt.1.) then
812c     print*,'WARN!!! ig=',ig,'  l=',l
813c    s         ,'   pdoadj=',pdoadj(ig,l)
814c           endif
815c        enddo
816c      enddo
817
818      print*,'14 OK convect8'
819c------------------------------------------------------------------
820c   Calculs pour les sorties
821c------------------------------------------------------------------
822
823      if(sorties) then
824      do l=1,nlay
825         do ig=1,ngrid
826            zla(ig,l)=(1.-fracd(ig,l))*zmax(ig)
827            zld(ig,l)=fracd(ig,l)*zmax(ig)
828            if(1.-fracd(ig,l).gt.1.e-10)
829     s      zwa(ig,l)=wd(ig,l)*fracd(ig,l)/(1.-fracd(ig,l))
830         enddo
831      enddo
832
833cdeja fait
834c      do l=1,nlay
835c         do ig=1,ngrid
836c            detr(ig,l)=fm(ig,l)+entr(ig,l)-fm(ig,l+1)
837c            if (detr(ig,l).lt.0.) then
838c                entr(ig,l)=entr(ig,l)-detr(ig,l)
839c                detr(ig,l)=0.
840c     print*,'WARNING !!! detrainement negatif ',ig,l
841c            endif
842c         enddo
843c      enddo
844
845c     print*,'15 OK convect8'
846
847      isplit=isplit+1
848
849
850c #define und
851        goto 123
852#ifdef und
853      CALL writeg1d(1,nlay,wd,'wd      ','wd      ')
854      CALL writeg1d(1,nlay,zwa,'wa      ','wa      ')
855      CALL writeg1d(1,nlay,fracd,'fracd      ','fracd      ')
856      CALL writeg1d(1,nlay,fraca,'fraca      ','fraca      ')
857      CALL writeg1d(1,nlay,wa_moy,'wam         ','wam         ')
858      CALL writeg1d(1,nlay,zla,'la      ','la      ')
859      CALL writeg1d(1,nlay,zld,'ld      ','ld      ')
860      CALL writeg1d(1,nlay,pt,'pt      ','pt      ')
861      CALL writeg1d(1,nlay,zh,'zh      ','zh      ')
862      CALL writeg1d(1,nlay,zha,'zha      ','zha      ')
863      CALL writeg1d(1,nlay,zu,'zu      ','zu      ')
864      CALL writeg1d(1,nlay,zv,'zv      ','zv      ')
865      CALL writeg1d(1,nlay,zo,'zo      ','zo      ')
866      CALL writeg1d(1,nlay,wh,'wh      ','wh      ')
867      CALL writeg1d(1,nlay,wu,'wu      ','wu      ')
868      CALL writeg1d(1,nlay,wv,'wv      ','wv      ')
869      CALL writeg1d(1,nlay,wo,'w15uo     ','wXo     ')
870      CALL writeg1d(1,nlay,zdhadj,'zdhadj      ','zdhadj      ')
871      CALL writeg1d(1,nlay,pduadj,'pduadj      ','pduadj      ')
872      CALL writeg1d(1,nlay,pdvadj,'pdvadj      ','pdvadj      ')
873      CALL writeg1d(1,nlay,pdoadj,'pdoadj      ','pdoadj      ')
874      CALL writeg1d(1,nlay,entr  ,'entr        ','entr        ')
875      CALL writeg1d(1,nlay,detr  ,'detr        ','detr        ')
876      CALL writeg1d(1,nlay,fm    ,'fm          ','fm          ')
877
878      CALL writeg1d(1,nlay,pdtadj,'pdtadj    ','pdtadj    ')
879      CALL writeg1d(1,nlay,pplay,'pplay     ','pplay     ')
880      CALL writeg1d(1,nlay,pplev,'pplev     ','pplev     ')
881
882c   recalcul des flux en diagnostique...
883c     print*,'PAS DE TEMPS ',ptimestep
884       call dt2F(pplev,pplay,pt,pdtadj,wh)
885      CALL writeg1d(1,nlay,wh,'wh2     ','wh2     ')
886#endif
887123   continue
888c#define troisD
889#ifdef troisD
890c       if (sorties) then
891      print*,'Debut des wrgradsfi'
892
893c      print*,'16 OK convect8'
894         call wrgradsfi(1,nlay,wd,'wd        ','wd        ')
895         call wrgradsfi(1,nlay,zwa,'wa        ','wa        ')
896         call wrgradsfi(1,nlay,fracd,'fracd     ','fracd     ')
897         call wrgradsfi(1,nlay,fraca,'fraca     ','fraca     ')
898         call wrgradsfi(1,nlay,xxx,'xxx       ','xxx       ')
899         call wrgradsfi(1,nlay,wa_moy,'wam       ','wam       ')
900c      print*,'WA6 ',wa_moy
901         call wrgradsfi(1,nlay,zla,'la        ','la        ')
902         call wrgradsfi(1,nlay,zld,'ld        ','ld        ')
903         call wrgradsfi(1,nlay,pt,'pt        ','pt        ')
904         call wrgradsfi(1,nlay,zh,'zh        ','zh        ')
905         call wrgradsfi(1,nlay,zha,'zha        ','zha        ')
906         call wrgradsfi(1,nlay,zua,'zua        ','zua        ')
907         call wrgradsfi(1,nlay,zva,'zva        ','zva        ')
908         call wrgradsfi(1,nlay,zu,'zu        ','zu        ')
909         call wrgradsfi(1,nlay,zv,'zv        ','zv        ')
910         call wrgradsfi(1,nlay,zo,'zo        ','zo        ')
911         call wrgradsfi(1,nlay,wh,'wh        ','wh        ')
912         call wrgradsfi(1,nlay,wu,'wu        ','wu        ')
913         call wrgradsfi(1,nlay,wv,'wv        ','wv        ')
914         call wrgradsfi(1,nlay,wo,'wo        ','wo        ')
915         call wrgradsfi(1,1,zmax,'zmax      ','zmax      ')
916         call wrgradsfi(1,nlay,zdhadj,'zdhadj    ','zdhadj    ')
917         call wrgradsfi(1,nlay,pduadj,'pduadj    ','pduadj    ')
918         call wrgradsfi(1,nlay,pdvadj,'pdvadj    ','pdvadj    ')
919         call wrgradsfi(1,nlay,pdoadj,'pdoadj    ','pdoadj    ')
920         call wrgradsfi(1,nlay,entr,'entr      ','entr      ')
921         call wrgradsfi(1,nlay,detr,'detr      ','detr      ')
922         call wrgradsfi(1,nlay,fm,'fm        ','fm        ')
923         call wrgradsfi(1,nlay,fmc,'fmc       ','fmc       ')
924         call wrgradsfi(1,nlay,zw2,'zw2       ','zw2       ')
925         call wrgradsfi(1,nlay,ztva,'ztva      ','ztva      ')
926         call wrgradsfi(1,nlay,ztv,'ztv       ','ztv       ')
927
928         call wrgradsfi(1,nlay,zo,'zo        ','zo        ')
929         call wrgradsfi(1,nlay,larg_cons,'Lc        ','Lc        ')
930         call wrgradsfi(1,nlay,larg_detr,'Ldetr     ','Ldetr     ')
931
932cCR:nouveaux diagnostiques
933      call wrgradsfi(1,nlay,entr_star  ,'entr_star   ','entr_star   ')     
934      call wrgradsfi(1,nlay,f_star    ,'f_star   ','f_star   ')
935      call wrgradsfi(1,1,zmax,'zmax      ','zmax      ')
936      call wrgradsfi(1,1,zmix,'zmix      ','zmix      ')
937      zsortie1d(:)=lmax(:)
938      call wrgradsfi(1,1,zsortie1d,'lmax      ','lmax      ')
939      call wrgradsfi(1,1,wmax,'wmax      ','wmax      ')
940      zsortie1d(:)=lmix(:)
941      call wrgradsfi(1,1,zsortie1d,'lmix      ','lmix      ')
942      zsortie1d(:)=lentr(:)
943      call wrgradsfi(1,1,zsortie1d,'lentr      ','lentr     ')
944
945c      print*,'17 OK convect8'
946
947         do k=1,klev/10
948            write(str2,'(i2.2)') k
949            str10='wa'//str2
950            do l=1,nlay
951               do ig=1,ngrid
952                  zsortie(ig,l)=wa(ig,k,l)
953               enddo
954            enddo   
955            CALL wrgradsfi(1,nlay,zsortie,str10,str10)
956            do l=1,nlay
957               do ig=1,ngrid
958                  zsortie(ig,l)=larg_part(ig,k,l)
959               enddo
960            enddo
961            str10='la'//str2
962            CALL wrgradsfi(1,nlay,zsortie,str10,str10)
963         enddo
964
965
966c     print*,'18 OK convect8'
967c      endif
968      print*,'Fin des wrgradsfi'
969#endif
970
971      endif
972
973c     if(wa_moy(1,4).gt.1.e-10) stop
974
975      print*,'19 OK convect8'
976      return
977      end
978
979      subroutine dqthermcell(ngrid,nlay,ptimestep,fm,entr,masse
980     .    ,q,dq,qa)
981      implicit none
982
983c=======================================================================
984c
985c   Calcul du transport verticale dans la couche limite en presence
986c   de "thermiques" explicitement representes
987c   calcul du dq/dt une fois qu'on connait les ascendances
988c
989c=======================================================================
990
991#include "dimensions.h"
992#include "dimphy.h"
993
994      integer ngrid,nlay
995
996      real ptimestep
997      real masse(ngrid,nlay),fm(ngrid,nlay+1)
998      real entr(ngrid,nlay)
999      real q(ngrid,nlay)
1000      real dq(ngrid,nlay)
1001
1002      real qa(klon,klev),detr(klon,klev),wqd(klon,klev+1)
1003
1004      integer ig,k
1005
1006c   calcul du detrainement
1007
1008      do k=1,nlay
1009         do ig=1,ngrid
1010            detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
1011         enddo
1012      enddo
1013
1014c   calcul de la valeur dans les ascendances
1015      do ig=1,ngrid
1016         qa(ig,1)=q(ig,1)
1017      enddo
1018
1019      do k=2,nlay
1020         do ig=1,ngrid
1021            if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.
1022     s         1.e-5*masse(ig,k)) then
1023               qa(ig,k)=(fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k))
1024     s         /(fm(ig,k+1)+detr(ig,k))
1025            else
1026               qa(ig,k)=q(ig,k)
1027            endif
1028         enddo
1029      enddo
1030
1031      do k=2,nlay
1032         do ig=1,ngrid
1033c             wqd(ig,k)=fm(ig,k)*0.5*(q(ig,k-1)+q(ig,k))
1034            wqd(ig,k)=fm(ig,k)*q(ig,k)
1035         enddo
1036      enddo
1037      do ig=1,ngrid
1038         wqd(ig,1)=0.
1039         wqd(ig,nlay+1)=0.
1040      enddo
1041
1042      do k=1,nlay
1043         do ig=1,ngrid
1044            dq(ig,k)=(detr(ig,k)*qa(ig,k)-entr(ig,k)*q(ig,k)
1045     s               -wqd(ig,k)+wqd(ig,k+1))
1046     s               /masse(ig,k)
1047         enddo
1048      enddo
1049
1050      return
1051      end
1052      subroutine dvthermcell(ngrid,nlay,ptimestep,fm,entr,masse
1053     .    ,fraca,larga
1054     .    ,u,v,du,dv,ua,va)
1055      implicit none
1056
1057c=======================================================================
1058c
1059c   Calcul du transport verticale dans la couche limite en presence
1060c   de "thermiques" explicitement representes
1061c   calcul du dq/dt une fois qu'on connait les ascendances
1062c
1063c=======================================================================
1064
1065#include "dimensions.h"
1066#include "dimphy.h"
1067
1068      integer ngrid,nlay
1069
1070      real ptimestep
1071      real masse(ngrid,nlay),fm(ngrid,nlay+1)
1072      real fraca(ngrid,nlay+1)
1073      real larga(ngrid)
1074      real entr(ngrid,nlay)
1075      real u(ngrid,nlay)
1076      real ua(ngrid,nlay)
1077      real du(ngrid,nlay)
1078      real v(ngrid,nlay)
1079      real va(ngrid,nlay)
1080      real dv(ngrid,nlay)
1081
1082      real qa(klon,klev),detr(klon,klev)
1083      real wvd(klon,klev+1),wud(klon,klev+1)
1084      real gamma0,gamma(klon,klev+1)
1085      real dua,dva
1086      integer iter
1087
1088      integer ig,k
1089
1090c   calcul du detrainement
1091
1092      do k=1,nlay
1093         do ig=1,ngrid
1094            detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
1095         enddo
1096      enddo
1097
1098c   calcul de la valeur dans les ascendances
1099      do ig=1,ngrid
1100         ua(ig,1)=u(ig,1)
1101         va(ig,1)=v(ig,1)
1102      enddo
1103
1104      do k=2,nlay
1105         do ig=1,ngrid
1106            if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.
1107     s         1.e-5*masse(ig,k)) then
1108c   On itère sur la valeur du coeff de freinage.
1109c              gamma0=rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k))
1110               gamma0=masse(ig,k)
1111     s         *sqrt( 0.5*(fraca(ig,k+1)+fraca(ig,k)) )
1112     s         *0.5/larga(ig)
1113c              gamma0=0.
1114c   la première fois on multiplie le coefficient de freinage
1115c   par le module du vent dans la couche en dessous.
1116               dua=ua(ig,k-1)-u(ig,k-1)
1117               dva=va(ig,k-1)-v(ig,k-1)
1118               do iter=1,5
1119                  gamma(ig,k)=gamma0*sqrt(dua**2+dva**2)
1120                  ua(ig,k)=(fm(ig,k)*ua(ig,k-1)
1121     s               +(entr(ig,k)+gamma(ig,k))*u(ig,k))
1122     s               /(fm(ig,k+1)+detr(ig,k)+gamma(ig,k))
1123                  va(ig,k)=(fm(ig,k)*va(ig,k-1)
1124     s               +(entr(ig,k)+gamma(ig,k))*v(ig,k))
1125     s               /(fm(ig,k+1)+detr(ig,k)+gamma(ig,k))
1126c                 print*,k,ua(ig,k),va(ig,k),u(ig,k),v(ig,k),dua,dva
1127                  dua=ua(ig,k)-u(ig,k)
1128                  dva=va(ig,k)-v(ig,k)
1129               enddo
1130            else
1131               ua(ig,k)=u(ig,k)
1132               va(ig,k)=v(ig,k)
1133               gamma(ig,k)=0.
1134            endif
1135         enddo
1136      enddo
1137
1138      do k=2,nlay
1139         do ig=1,ngrid
1140            wud(ig,k)=fm(ig,k)*u(ig,k)
1141            wvd(ig,k)=fm(ig,k)*v(ig,k)
1142         enddo
1143      enddo
1144      do ig=1,ngrid
1145         wud(ig,1)=0.
1146         wud(ig,nlay+1)=0.
1147         wvd(ig,1)=0.
1148         wvd(ig,nlay+1)=0.
1149      enddo
1150
1151      do k=1,nlay
1152         do ig=1,ngrid
1153            du(ig,k)=((detr(ig,k)+gamma(ig,k))*ua(ig,k)
1154     s               -(entr(ig,k)+gamma(ig,k))*u(ig,k)
1155     s               -wud(ig,k)+wud(ig,k+1))
1156     s               /masse(ig,k)
1157            dv(ig,k)=((detr(ig,k)+gamma(ig,k))*va(ig,k)
1158     s               -(entr(ig,k)+gamma(ig,k))*v(ig,k)
1159     s               -wvd(ig,k)+wvd(ig,k+1))
1160     s               /masse(ig,k)
1161         enddo
1162      enddo
1163
1164      return
1165      end
1166      subroutine dqthermcell2(ngrid,nlay,ptimestep,fm,entr,masse,frac
1167     .    ,q,dq,qa)
1168      implicit none
1169
1170c=======================================================================
1171c
1172c   Calcul du transport verticale dans la couche limite en presence
1173c   de "thermiques" explicitement representes
1174c   calcul du dq/dt une fois qu'on connait les ascendances
1175c
1176c=======================================================================
1177
1178#include "dimensions.h"
1179#include "dimphy.h"
1180
1181      integer ngrid,nlay
1182
1183      real ptimestep
1184      real masse(ngrid,nlay),fm(ngrid,nlay+1)
1185      real entr(ngrid,nlay),frac(ngrid,nlay)
1186      real q(ngrid,nlay)
1187      real dq(ngrid,nlay)
1188
1189      real qa(klon,klev),detr(klon,klev),wqd(klon,klev+1)
1190      real qe(klon,klev),zf,zf2
1191
1192      integer ig,k
1193
1194c   calcul du detrainement
1195
1196      do k=1,nlay
1197         do ig=1,ngrid
1198            detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
1199         enddo
1200      enddo
1201
1202c   calcul de la valeur dans les ascendances
1203      do ig=1,ngrid
1204         qa(ig,1)=q(ig,1)
1205         qe(ig,1)=q(ig,1)
1206      enddo
1207
1208      do k=2,nlay
1209         do ig=1,ngrid
1210            if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.
1211     s         1.e-5*masse(ig,k)) then
1212               zf=0.5*(frac(ig,k)+frac(ig,k+1))
1213               zf2=1./(1.-zf)
1214               qa(ig,k)=(fm(ig,k)*qa(ig,k-1)+zf2*entr(ig,k)*q(ig,k))
1215     s         /(fm(ig,k+1)+detr(ig,k)+entr(ig,k)*zf*zf2)
1216               qe(ig,k)=(q(ig,k)-zf*qa(ig,k))*zf2
1217            else
1218               qa(ig,k)=q(ig,k)
1219               qe(ig,k)=q(ig,k)
1220            endif
1221         enddo
1222      enddo
1223
1224      do k=2,nlay
1225         do ig=1,ngrid
1226c             wqd(ig,k)=fm(ig,k)*0.5*(q(ig,k-1)+q(ig,k))
1227            wqd(ig,k)=fm(ig,k)*qe(ig,k)
1228         enddo
1229      enddo
1230      do ig=1,ngrid
1231         wqd(ig,1)=0.
1232         wqd(ig,nlay+1)=0.
1233      enddo
1234
1235      do k=1,nlay
1236         do ig=1,ngrid
1237            dq(ig,k)=(detr(ig,k)*qa(ig,k)-entr(ig,k)*qe(ig,k)
1238     s               -wqd(ig,k)+wqd(ig,k+1))
1239     s               /masse(ig,k)
1240         enddo
1241      enddo
1242
1243      return
1244      end
1245      subroutine dvthermcell2(ngrid,nlay,ptimestep,fm,entr,masse
1246     .    ,fraca,larga
1247     .    ,u,v,du,dv,ua,va)
1248      implicit none
1249
1250c=======================================================================
1251c
1252c   Calcul du transport verticale dans la couche limite en presence
1253c   de "thermiques" explicitement representes
1254c   calcul du dq/dt une fois qu'on connait les ascendances
1255c
1256c=======================================================================
1257
1258#include "dimensions.h"
1259#include "dimphy.h"
1260
1261      integer ngrid,nlay
1262
1263      real ptimestep
1264      real masse(ngrid,nlay),fm(ngrid,nlay+1)
1265      real fraca(ngrid,nlay+1)
1266      real larga(ngrid)
1267      real entr(ngrid,nlay)
1268      real u(ngrid,nlay)
1269      real ua(ngrid,nlay)
1270      real du(ngrid,nlay)
1271      real v(ngrid,nlay)
1272      real va(ngrid,nlay)
1273      real dv(ngrid,nlay)
1274
1275      real qa(klon,klev),detr(klon,klev),zf,zf2
1276      real wvd(klon,klev+1),wud(klon,klev+1)
1277      real gamma0,gamma(klon,klev+1)
1278      real ue(klon,klev),ve(klon,klev)
1279      real dua,dva
1280      integer iter
1281
1282      integer ig,k
1283
1284c   calcul du detrainement
1285
1286      do k=1,nlay
1287         do ig=1,ngrid
1288            detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
1289         enddo
1290      enddo
1291
1292c   calcul de la valeur dans les ascendances
1293      do ig=1,ngrid
1294         ua(ig,1)=u(ig,1)
1295         va(ig,1)=v(ig,1)
1296         ue(ig,1)=u(ig,1)
1297         ve(ig,1)=v(ig,1)
1298      enddo
1299
1300      do k=2,nlay
1301         do ig=1,ngrid
1302            if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.
1303     s         1.e-5*masse(ig,k)) then
1304c   On itère sur la valeur du coeff de freinage.
1305c              gamma0=rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k))
1306               gamma0=masse(ig,k)
1307     s         *sqrt( 0.5*(fraca(ig,k+1)+fraca(ig,k)) )
1308     s         *0.5/larga(ig)
1309     s         *1.
1310c    s         *0.5
1311c              gamma0=0.
1312               zf=0.5*(fraca(ig,k)+fraca(ig,k+1))
1313               zf=0.
1314               zf2=1./(1.-zf)
1315c   la première fois on multiplie le coefficient de freinage
1316c   par le module du vent dans la couche en dessous.
1317               dua=ua(ig,k-1)-u(ig,k-1)
1318               dva=va(ig,k-1)-v(ig,k-1)
1319               do iter=1,5
1320c   On choisit une relaxation lineaire.
1321                  gamma(ig,k)=gamma0
1322c   On choisit une relaxation quadratique.
1323                  gamma(ig,k)=gamma0*sqrt(dua**2+dva**2)
1324                  ua(ig,k)=(fm(ig,k)*ua(ig,k-1)
1325     s               +(zf2*entr(ig,k)+gamma(ig,k))*u(ig,k))
1326     s               /(fm(ig,k+1)+detr(ig,k)+entr(ig,k)*zf*zf2
1327     s                 +gamma(ig,k))
1328                  va(ig,k)=(fm(ig,k)*va(ig,k-1)
1329     s               +(zf2*entr(ig,k)+gamma(ig,k))*v(ig,k))
1330     s               /(fm(ig,k+1)+detr(ig,k)+entr(ig,k)*zf*zf2
1331     s                 +gamma(ig,k))
1332c                 print*,k,ua(ig,k),va(ig,k),u(ig,k),v(ig,k),dua,dva
1333                  dua=ua(ig,k)-u(ig,k)
1334                  dva=va(ig,k)-v(ig,k)
1335                  ue(ig,k)=(u(ig,k)-zf*ua(ig,k))*zf2
1336                  ve(ig,k)=(v(ig,k)-zf*va(ig,k))*zf2
1337               enddo
1338            else
1339               ua(ig,k)=u(ig,k)
1340               va(ig,k)=v(ig,k)
1341               ue(ig,k)=u(ig,k)
1342               ve(ig,k)=v(ig,k)
1343               gamma(ig,k)=0.
1344            endif
1345         enddo
1346      enddo
1347
1348      do k=2,nlay
1349         do ig=1,ngrid
1350            wud(ig,k)=fm(ig,k)*ue(ig,k)
1351            wvd(ig,k)=fm(ig,k)*ve(ig,k)
1352         enddo
1353      enddo
1354      do ig=1,ngrid
1355         wud(ig,1)=0.
1356         wud(ig,nlay+1)=0.
1357         wvd(ig,1)=0.
1358         wvd(ig,nlay+1)=0.
1359      enddo
1360
1361      do k=1,nlay
1362         do ig=1,ngrid
1363            du(ig,k)=((detr(ig,k)+gamma(ig,k))*ua(ig,k)
1364     s               -(entr(ig,k)+gamma(ig,k))*ue(ig,k)
1365     s               -wud(ig,k)+wud(ig,k+1))
1366     s               /masse(ig,k)
1367            dv(ig,k)=((detr(ig,k)+gamma(ig,k))*va(ig,k)
1368     s               -(entr(ig,k)+gamma(ig,k))*ve(ig,k)
1369     s               -wvd(ig,k)+wvd(ig,k+1))
1370     s               /masse(ig,k)
1371         enddo
1372      enddo
1373
1374      return
1375      end
Note: See TracBrowser for help on using the repository browser.