source: LMDZ4/branches/V3_test/libf/phylmd/thermcell.F @ 2116

Last change on this file since 2116 was 735, checked in by Laurent Fairhead, 18 years ago

Modifications version parallele
YM/LF

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