source: LMDZ5/trunk/libf/phylmd/thermcell.F90 @ 2415

Last change on this file since 2415 was 2346, checked in by Ehouarn Millour, 9 years ago

Physics/dynamics separation:

  • remove all references to dimensions.h from physics. nbp_lon (==iim) , nbp_lat (==jjm+1) and nbp_lev (==llm) from mod_grid_phy_lmdz should be used instead.
  • added module regular_lonlat_mod in phy_common to store information about the global (lon-lat) grid cell boundaries and centers.

EM
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 32.5 KB
Line 
1
2! $Id: thermcell.F90 2346 2015-08-21 15:13:46Z crio $
3
4SUBROUTINE calcul_sec(ngrid, nlay, ptimestep, pplay, pplev, pphi, zlev, pu, &
5    pv, pt, po, zmax, wmax, zw2, lmix & ! s
6                                        ! ,pu_therm,pv_therm
7    , r_aspect, l_mix, w2di, tho)
8
9  USE dimphy
10  IMPLICIT NONE
11
12  ! =======================================================================
13
14  ! Calcul du transport verticale dans la couche limite en presence
15  ! de "thermiques" explicitement representes
16
17  ! Réécriture à partir d'un listing papier à Habas, le 14/02/00
18
19  ! le thermique est supposé homogène et dissipé par mélange avec
20  ! son environnement. la longueur l_mix contrôle l'efficacité du
21  ! mélange
22
23  ! Le calcul du transport des différentes espèces se fait en prenant
24  ! en compte:
25  ! 1. un flux de masse montant
26  ! 2. un flux de masse descendant
27  ! 3. un entrainement
28  ! 4. un detrainement
29
30  ! =======================================================================
31
32  ! -----------------------------------------------------------------------
33  ! declarations:
34  ! -------------
35
36  include "YOMCST.h"
37
38  ! arguments:
39  ! ----------
40
41  INTEGER ngrid, nlay, w2di
42  REAL 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  !$OMP THREADPRIVATE(idetr)
55  ! local:
56  ! ------
57
58  INTEGER ig, k, l, lmaxa(klon), lmix(klon)
59  REAL zsortie1d(klon)
60  ! 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)
64  ! RC
65  REAL zmax(klon), zw, zw2(klon, klev+1), ztva(klon, klev)
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
92  INTEGER isplit, nsplit, ialt
93  PARAMETER (nsplit=10)
94  DATA isplit/0/
95  SAVE isplit
96  !$OMP THREADPRIVATE(isplit)
97
98  LOGICAL sorties
99  REAL rho(klon, klev), rhobarz(klon, klev+1), masse(klon, klev)
100  REAL zpspsk(klon, klev)
101
102  ! real wmax(klon,klev),wmaxa(klon)
103  REAL wmax(klon), wmaxa(klon)
104  REAL wa(klon, klev, klev+1)
105  REAL wd(klon, klev+1)
106  REAL larg_part(klon, klev, klev+1)
107  REAL fracd(klon, klev+1)
108  REAL xxx(klon, klev+1)
109  REAL larg_cons(klon, klev+1)
110  REAL larg_detr(klon, klev+1)
111  REAL fm0(klon, klev+1), entr0(klon, klev), detr(klon, klev)
112  REAL pu_therm(klon, klev), pv_therm(klon, klev)
113  REAL fm(klon, klev+1), entr(klon, klev)
114  REAL fmc(klon, klev+1)
115
116  ! CR:nouvelles variables
117  REAL f_star(klon, klev+1), entr_star(klon, klev)
118  REAL entr_star_tot(klon), entr_star2(klon)
119  REAL zalim(klon)
120  INTEGER lalim(klon)
121  REAL norme(klon)
122  REAL f(klon), f0(klon)
123  REAL zlevinter(klon)
124  LOGICAL therm
125  LOGICAL first
126  DATA first/.FALSE./
127  SAVE first
128  !$OMP THREADPRIVATE(first)
129  ! RC
130
131  CHARACTER *2 str2
132  CHARACTER *10 str10
133
134  CHARACTER (LEN=20) :: modname = 'calcul_sec'
135  CHARACTER (LEN=80) :: abort_message
136
137
138  ! LOGICAL vtest(klon),down
139
140  EXTERNAL scopy
141
142  INTEGER ncorrec
143  SAVE ncorrec
144  DATA ncorrec/0/
145  !$OMP THREADPRIVATE(ncorrec)
146
147
148  ! -----------------------------------------------------------------------
149  ! initialisation:
150  ! ---------------
151
152  sorties = .TRUE.
153  IF (ngrid/=klon) THEN
154    PRINT *
155    PRINT *, 'STOP dans convadj'
156    PRINT *, 'ngrid    =', ngrid
157    PRINT *, 'klon  =', klon
158  END IF
159
160  ! -----------------------------------------------------------------------
161  ! incrementation eventuelle de tendances precedentes:
162  ! ---------------------------------------------------
163
164  ! print*,'0 OK convect8'
165
166  DO l = 1, nlay
167    DO ig = 1, ngrid
168      zpspsk(ig, l) = (pplay(ig,l)/pplev(ig,1))**rkappa
169      zh(ig, l) = pt(ig, l)/zpspsk(ig, l)
170      zu(ig, l) = pu(ig, l)
171      zv(ig, l) = pv(ig, l)
172      zo(ig, l) = po(ig, l)
173      ztv(ig, l) = zh(ig, l)*(1.+0.61*zo(ig,l))
174    END DO
175  END DO
176
177  ! print*,'1 OK convect8'
178  ! --------------------
179
180
181  ! + + + + + + + + + + +
182
183
184  ! wa, fraca, wd, fracd --------------------   zlev(2), rhobarz
185  ! wh,wt,wo ...
186
187  ! + + + + + + + + + + +  zh,zu,zv,zo,rho
188
189
190  ! --------------------   zlev(1)
191  ! \\\\\\\\\\\\\\\\\\\\
192
193
194
195  ! -----------------------------------------------------------------------
196  ! Calcul des altitudes des couches
197  ! -----------------------------------------------------------------------
198
199  DO l = 2, nlay
200    DO ig = 1, ngrid
201      zlev(ig, l) = 0.5*(pphi(ig,l)+pphi(ig,l-1))/rg
202    END DO
203  END DO
204  DO ig = 1, ngrid
205    zlev(ig, 1) = 0.
206    zlev(ig, nlay+1) = (2.*pphi(ig,klev)-pphi(ig,klev-1))/rg
207  END DO
208  DO l = 1, nlay
209    DO ig = 1, ngrid
210      zlay(ig, l) = pphi(ig, l)/rg
211    END DO
212  END DO
213
214  ! print*,'2 OK convect8'
215  ! -----------------------------------------------------------------------
216  ! Calcul des densites
217  ! -----------------------------------------------------------------------
218
219  DO l = 1, nlay
220    DO ig = 1, ngrid
221      rho(ig, l) = pplay(ig, l)/(zpspsk(ig,l)*rd*zh(ig,l))
222    END DO
223  END DO
224
225  DO l = 2, nlay
226    DO ig = 1, ngrid
227      rhobarz(ig, l) = 0.5*(rho(ig,l)+rho(ig,l-1))
228    END DO
229  END DO
230
231  DO k = 1, nlay
232    DO l = 1, nlay + 1
233      DO ig = 1, ngrid
234        wa(ig, k, l) = 0.
235      END DO
236    END DO
237  END DO
238
239  ! print*,'3 OK convect8'
240  ! ------------------------------------------------------------------
241  ! Calcul de w2, quarre de w a partir de la cape
242  ! a partir de w2, on calcule wa, vitesse de l'ascendance
243
244  ! ATTENTION: Dans cette version, pour cause d'economie de memoire,
245  ! w2 est stoke dans wa
246
247  ! ATTENTION: dans convect8, on n'utilise le calcule des wa
248  ! independants par couches que pour calculer l'entrainement
249  ! a la base et la hauteur max de l'ascendance.
250
251  ! Indicages:
252  ! l'ascendance provenant du niveau k traverse l'interface l avec
253  ! une vitesse wa(k,l).
254
255  ! --------------------
256
257  ! + + + + + + + + + +
258
259  ! wa(k,l)   ----       --------------------    l
260  ! /\
261  ! /||\       + + + + + + + + + +
262  ! ||
263  ! ||        --------------------
264  ! ||
265  ! ||        + + + + + + + + + +
266  ! ||
267  ! ||        --------------------
268  ! ||__
269  ! |___      + + + + + + + + + +     k
270
271  ! --------------------
272
273
274
275  ! ------------------------------------------------------------------
276
277  ! CR: ponderation entrainement des couches instables
278  ! def des entr_star tels que entr=f*entr_star
279  DO l = 1, klev
280    DO ig = 1, ngrid
281      entr_star(ig, l) = 0.
282    END DO
283  END DO
284  ! determination de la longueur de la couche d entrainement
285  DO ig = 1, ngrid
286    lentr(ig) = 1
287  END DO
288
289  ! on ne considere que les premieres couches instables
290  therm = .FALSE.
291  DO k = nlay - 2, 1, -1
292    DO ig = 1, ngrid
293      IF (ztv(ig,k)>ztv(ig,k+1) .AND. ztv(ig,k+1)<=ztv(ig,k+2)) THEN
294        lentr(ig) = k + 1
295        therm = .TRUE.
296      END IF
297    END DO
298  END DO
299  ! limitation de la valeur du lentr
300  ! do ig=1,ngrid
301  ! lentr(ig)=min(5,lentr(ig))
302  ! enddo
303  ! determination du lmin: couche d ou provient le thermique
304  DO ig = 1, ngrid
305    lmin(ig) = 1
306  END DO
307  DO ig = 1, ngrid
308    DO l = nlay, 2, -1
309      IF (ztv(ig,l-1)>ztv(ig,l)) THEN
310        lmin(ig) = l - 1
311      END IF
312    END DO
313  END DO
314  ! initialisations
315  DO ig = 1, ngrid
316    zalim(ig) = 0.
317    norme(ig) = 0.
318    lalim(ig) = 1
319  END DO
320  DO k = 1, klev - 1
321    DO ig = 1, ngrid
322      zalim(ig) = zalim(ig) + zlev(ig, k)*max(0., (ztv(ig,k)-ztv(ig, &
323        k+1))/(zlev(ig,k+1)-zlev(ig,k)))
324      ! s         *(zlev(ig,k+1)-zlev(ig,k))
325      norme(ig) = norme(ig) + max(0., (ztv(ig,k)-ztv(ig,k+1))/(zlev(ig, &
326        k+1)-zlev(ig,k)))
327      ! s          *(zlev(ig,k+1)-zlev(ig,k))
328    END DO
329  END DO
330  DO ig = 1, ngrid
331    IF (norme(ig)>1.E-10) THEN
332      zalim(ig) = max(10.*zalim(ig)/norme(ig), zlev(ig,2))
333      ! zalim(ig)=min(zalim(ig),zlev(ig,lentr(ig)))
334    END IF
335  END DO
336  ! détermination du lalim correspondant
337  DO k = 1, klev - 1
338    DO ig = 1, ngrid
339      IF ((zalim(ig)>zlev(ig,k)) .AND. (zalim(ig)<=zlev(ig,k+1))) THEN
340        lalim(ig) = k
341      END IF
342    END DO
343  END DO
344
345  ! definition de l'entrainement des couches
346  DO l = 1, klev - 1
347    DO ig = 1, ngrid
348      IF (ztv(ig,l)>ztv(ig,l+1) .AND. l>=lmin(ig) .AND. l<lentr(ig)) THEN
349        entr_star(ig, l) = max((ztv(ig,l)-ztv(ig,l+1)), 0.) & ! s
350                                                              ! *(zlev(ig,l+1)-zlev(ig,l))
351          *sqrt(zlev(ig,l+1))
352        ! autre def
353        ! entr_star(ig,l)=zlev(ig,l+1)*(1.-(zlev(ig,l+1)
354        ! s                         /zlev(ig,lentr(ig)+2)))**(3./2.)
355      END IF
356    END DO
357  END DO
358  ! nouveau test
359  ! if (therm) then
360  DO l = 1, klev - 1
361    DO ig = 1, ngrid
362      IF (ztv(ig,l)>ztv(ig,l+1) .AND. l>=lmin(ig) .AND. l<=lalim(ig) .AND. &
363          zalim(ig)>1.E-10) THEN
364        ! if (l.le.lentr(ig)) then
365        ! entr_star(ig,l)=zlev(ig,l+1)*(1.-(zlev(ig,l+1)
366        ! s                         /zalim(ig)))**(3./2.)
367        ! write(10,*)zlev(ig,l),entr_star(ig,l)
368      END IF
369    END DO
370  END DO
371  ! endif
372  ! pas de thermique si couche 1 stable
373  DO ig = 1, ngrid
374    IF (lmin(ig)>5) THEN
375      DO l = 1, klev
376        entr_star(ig, l) = 0.
377      END DO
378    END IF
379  END DO
380  ! calcul de l entrainement total
381  DO ig = 1, ngrid
382    entr_star_tot(ig) = 0.
383  END DO
384  DO ig = 1, ngrid
385    DO k = 1, klev
386      entr_star_tot(ig) = entr_star_tot(ig) + entr_star(ig, k)
387    END DO
388  END DO
389  ! Calcul entrainement normalise
390  DO ig = 1, ngrid
391    IF (entr_star_tot(ig)>1.E-10) THEN
392      ! do l=1,lentr(ig)
393      DO l = 1, klev
394        ! def possibles pour entr_star: zdthetadz, dthetadz, zdtheta
395        entr_star(ig, l) = entr_star(ig, l)/entr_star_tot(ig)
396      END DO
397    END IF
398  END DO
399
400  ! print*,'fin calcul entr_star'
401  DO k = 1, klev
402    DO ig = 1, ngrid
403      ztva(ig, k) = ztv(ig, k)
404    END DO
405  END DO
406  ! RC
407  ! print*,'7 OK convect8'
408  DO k = 1, klev + 1
409    DO ig = 1, ngrid
410      zw2(ig, k) = 0.
411      fmc(ig, k) = 0.
412      ! CR
413      f_star(ig, k) = 0.
414      ! RC
415      larg_cons(ig, k) = 0.
416      larg_detr(ig, k) = 0.
417      wa_moy(ig, k) = 0.
418    END DO
419  END DO
420
421  ! print*,'8 OK convect8'
422  DO ig = 1, ngrid
423    linter(ig) = 1.
424    lmaxa(ig) = 1
425    lmix(ig) = 1
426    wmaxa(ig) = 0.
427  END DO
428
429  ! CR:
430  DO l = 1, nlay - 2
431    DO ig = 1, ngrid
432      IF (ztv(ig,l)>ztv(ig,l+1) .AND. entr_star(ig,l)>1.E-10 .AND. &
433          zw2(ig,l)<1E-10) THEN
434        f_star(ig, l+1) = entr_star(ig, l)
435        ! test:calcul de dteta
436        zw2(ig, l+1) = 2.*rg*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig, l+1)* &
437          (zlev(ig,l+1)-zlev(ig,l))*0.4*pphi(ig, l)/(pphi(ig,l+1)-pphi(ig,l))
438        larg_detr(ig, l) = 0.
439      ELSE IF ((zw2(ig,l)>=1E-10) .AND. (f_star(ig,l)+entr_star(ig, &
440          l)>1.E-10)) THEN
441        f_star(ig, l+1) = f_star(ig, l) + entr_star(ig, l)
442        ztva(ig, l) = (f_star(ig,l)*ztva(ig,l-1)+entr_star(ig,l)*ztv(ig,l))/ &
443          f_star(ig, l+1)
444        zw2(ig, l+1) = zw2(ig, l)*(f_star(ig,l)/f_star(ig,l+1))**2 + &
445          2.*rg*(ztva(ig,l)-ztv(ig,l))/ztv(ig, l)*(zlev(ig,l+1)-zlev(ig,l))
446      END IF
447      ! determination de zmax continu par interpolation lineaire
448      IF (zw2(ig,l+1)<0.) THEN
449        ! test
450        IF (abs(zw2(ig,l+1)-zw2(ig,l))<1E-10) THEN
451          ! print*,'pb linter'
452        END IF
453        linter(ig) = (l*(zw2(ig,l+1)-zw2(ig,l))-zw2(ig,l))/(zw2(ig,l+1)-zw2( &
454          ig,l))
455        zw2(ig, l+1) = 0.
456        lmaxa(ig) = l
457      ELSE
458        IF (zw2(ig,l+1)<0.) THEN
459          ! print*,'pb1 zw2<0'
460        END IF
461        wa_moy(ig, l+1) = sqrt(zw2(ig,l+1))
462      END IF
463      IF (wa_moy(ig,l+1)>wmaxa(ig)) THEN
464        ! lmix est le niveau de la couche ou w (wa_moy) est maximum
465        lmix(ig) = l + 1
466        wmaxa(ig) = wa_moy(ig, l+1)
467      END IF
468    END DO
469  END DO
470  ! print*,'fin calcul zw2'
471
472  ! Calcul de la couche correspondant a la hauteur du thermique
473  DO ig = 1, ngrid
474    lmax(ig) = lentr(ig)
475    ! lmax(ig)=lalim(ig)
476  END DO
477  DO ig = 1, ngrid
478    DO l = nlay, lentr(ig) + 1, -1
479      ! do l=nlay,lalim(ig)+1,-1
480      IF (zw2(ig,l)<=1.E-10) THEN
481        lmax(ig) = l - 1
482      END IF
483    END DO
484  END DO
485  ! pas de thermique si couche 1 stable
486  DO ig = 1, ngrid
487    IF (lmin(ig)>5) THEN
488      lmax(ig) = 1
489      lmin(ig) = 1
490      lentr(ig) = 1
491      lalim(ig) = 1
492    END IF
493  END DO
494
495  ! Determination de zw2 max
496  DO ig = 1, ngrid
497    wmax(ig) = 0.
498  END DO
499
500  DO l = 1, nlay
501    DO ig = 1, ngrid
502      IF (l<=lmax(ig)) THEN
503        IF (zw2(ig,l)<0.) THEN
504          ! print*,'pb2 zw2<0'
505        END IF
506        zw2(ig, l) = sqrt(zw2(ig,l))
507        wmax(ig) = max(wmax(ig), zw2(ig,l))
508      ELSE
509        zw2(ig, l) = 0.
510      END IF
511    END DO
512  END DO
513
514  ! Longueur caracteristique correspondant a la hauteur des thermiques.
515  DO ig = 1, ngrid
516    zmax(ig) = 0.
517    zlevinter(ig) = zlev(ig, 1)
518  END DO
519  DO ig = 1, ngrid
520    ! calcul de zlevinter
521    zlevinter(ig) = (zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*linter(ig) + &
522      zlev(ig, lmax(ig)) - lmax(ig)*(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))
523    zmax(ig) = max(zmax(ig), zlevinter(ig)-zlev(ig,lmin(ig)))
524  END DO
525  DO ig = 1, ngrid
526    ! write(8,*)zmax(ig),lmax(ig),lentr(ig),lmin(ig)
527  END DO
528  ! on stope après les calculs de zmax et wmax
529  RETURN
530
531  ! print*,'avant fermeture'
532  ! Fermeture,determination de f
533  ! Attention! entrainement normalisé ou pas?
534  DO ig = 1, ngrid
535    entr_star2(ig) = 0.
536  END DO
537  DO ig = 1, ngrid
538    IF (entr_star_tot(ig)<1.E-10) THEN
539      f(ig) = 0.
540    ELSE
541      DO k = lmin(ig), lentr(ig)
542        ! do k=lmin(ig),lalim(ig)
543        entr_star2(ig) = entr_star2(ig) + entr_star(ig, k)**2/(rho(ig,k)*( &
544          zlev(ig,k+1)-zlev(ig,k)))
545      END DO
546      ! Nouvelle fermeture
547      f(ig) = wmax(ig)/(max(500.,zmax(ig))*r_aspect*entr_star2(ig))
548      ! s            *entr_star_tot(ig)
549      ! test
550      ! if (first) then
551      f(ig) = f(ig) + (f0(ig)-f(ig))*exp(-ptimestep/zmax(ig)*wmax(ig))
552      ! endif
553    END IF
554    f0(ig) = f(ig)
555    ! first=.true.
556  END DO
557  ! print*,'apres fermeture'
558  ! on stoppe après la fermeture
559  RETURN
560  ! Calcul de l'entrainement
561  DO k = 1, klev
562    DO ig = 1, ngrid
563      entr(ig, k) = f(ig)*entr_star(ig, k)
564    END DO
565  END DO
566  ! on stoppe après le calcul de entr
567  ! RETURN
568  ! CR:test pour entrainer moins que la masse
569  ! do ig=1,ngrid
570  ! do l=1,lentr(ig)
571  ! if ((entr(ig,l)*ptimestep).gt.(0.9*masse(ig,l))) then
572  ! entr(ig,l+1)=entr(ig,l+1)+entr(ig,l)
573  ! s                       -0.9*masse(ig,l)/ptimestep
574  ! entr(ig,l)=0.9*masse(ig,l)/ptimestep
575  ! endif
576  ! enddo
577  ! enddo
578  ! CR: fin test
579  ! Calcul des flux
580  DO ig = 1, ngrid
581    DO l = 1, lmax(ig) - 1
582      fmc(ig, l+1) = fmc(ig, l) + entr(ig, l)
583    END DO
584  END DO
585
586  ! RC
587
588
589  ! print*,'9 OK convect8'
590  ! print*,'WA1 ',wa_moy
591
592  ! determination de l'indice du debut de la mixed layer ou w decroit
593
594  ! calcul de la largeur de chaque ascendance dans le cas conservatif.
595  ! dans ce cas simple, on suppose que la largeur de l'ascendance provenant
596  ! d'une couche est égale à la hauteur de la couche alimentante.
597  ! La vitesse maximale dans l'ascendance est aussi prise comme estimation
598  ! de la vitesse d'entrainement horizontal dans la couche alimentante.
599
600  DO l = 2, nlay
601    DO ig = 1, ngrid
602      IF (l<=lmaxa(ig)) THEN
603        zw = max(wa_moy(ig,l), 1.E-10)
604        larg_cons(ig, l) = zmax(ig)*r_aspect*fmc(ig, l)/(rhobarz(ig,l)*zw)
605      END IF
606    END DO
607  END DO
608
609  DO l = 2, nlay
610    DO ig = 1, ngrid
611      IF (l<=lmaxa(ig)) THEN
612        ! if (idetr.eq.0) then
613        ! cette option est finalement en dur.
614        IF ((l_mix*zlev(ig,l))<0.) THEN
615          ! print*,'pb l_mix*zlev<0'
616        END IF
617        ! CR: test: nouvelle def de lambda
618        ! larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
619        IF (zw2(ig,l)>1.E-10) THEN
620          larg_detr(ig, l) = sqrt((l_mix/zw2(ig,l))*zlev(ig,l))
621        ELSE
622          larg_detr(ig, l) = sqrt(l_mix*zlev(ig,l))
623        END IF
624        ! RC
625        ! else if (idetr.eq.1) then
626        ! larg_detr(ig,l)=larg_cons(ig,l)
627        ! s            *sqrt(l_mix*zlev(ig,l))/larg_cons(ig,lmix(ig))
628        ! else if (idetr.eq.2) then
629        ! larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
630        ! s            *sqrt(wa_moy(ig,l))
631        ! else if (idetr.eq.4) then
632        ! larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
633        ! s            *wa_moy(ig,l)
634        ! endif
635      END IF
636    END DO
637  END DO
638
639  ! print*,'10 OK convect8'
640  ! print*,'WA2 ',wa_moy
641  ! calcul de la fraction de la maille concernée par l'ascendance en tenant
642  ! compte de l'epluchage du thermique.
643
644  ! CR def de  zmix continu (profil parabolique des vitesses)
645  DO ig = 1, ngrid
646    IF (lmix(ig)>1.) THEN
647      ! test
648      IF (((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))*((zlev(ig,lmix(ig)))- &
649          (zlev(ig,lmix(ig)+1)))-(zw2(ig,lmix(ig))- &
650          zw2(ig,lmix(ig)+1))*((zlev(ig,lmix(ig)-1))- &
651          (zlev(ig,lmix(ig)))))>1E-10) THEN
652
653        zmix(ig) = ((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))*((zlev(ig,lmix(ig)) &
654          )**2-(zlev(ig,lmix(ig)+1))**2)-(zw2(ig,lmix(ig))-zw2(ig, &
655          lmix(ig)+1))*((zlev(ig,lmix(ig)-1))**2-(zlev(ig,lmix(ig)))**2))/ &
656          (2.*((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))*((zlev(ig,lmix(ig)))- &
657          (zlev(ig,lmix(ig)+1)))-(zw2(ig,lmix(ig))- &
658          zw2(ig,lmix(ig)+1))*((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))))
659      ELSE
660        zmix(ig) = zlev(ig, lmix(ig))
661        ! print*,'pb zmix'
662      END IF
663    ELSE
664      zmix(ig) = 0.
665    END IF
666    ! test
667    IF ((zmax(ig)-zmix(ig))<0.) THEN
668      zmix(ig) = 0.99*zmax(ig)
669      ! print*,'pb zmix>zmax'
670    END IF
671  END DO
672
673  ! calcul du nouveau lmix correspondant
674  DO ig = 1, ngrid
675    DO l = 1, klev
676      IF (zmix(ig)>=zlev(ig,l) .AND. zmix(ig)<zlev(ig,l+1)) THEN
677        lmix(ig) = l
678      END IF
679    END DO
680  END DO
681
682  DO l = 2, nlay
683    DO ig = 1, ngrid
684      IF (larg_cons(ig,l)>1.) THEN
685        ! print*,ig,l,lmix(ig),lmaxa(ig),larg_cons(ig,l),'  KKK'
686        fraca(ig, l) = (larg_cons(ig,l)-larg_detr(ig,l))/(r_aspect*zmax(ig))
687        ! test
688        fraca(ig, l) = max(fraca(ig,l), 0.)
689        fraca(ig, l) = min(fraca(ig,l), 0.5)
690        fracd(ig, l) = 1. - fraca(ig, l)
691        fracc(ig, l) = larg_cons(ig, l)/(r_aspect*zmax(ig))
692      ELSE
693        ! wa_moy(ig,l)=0.
694        fraca(ig, l) = 0.
695        fracc(ig, l) = 0.
696        fracd(ig, l) = 1.
697      END IF
698    END DO
699  END DO
700  ! CR: calcul de fracazmix
701  DO ig = 1, ngrid
702    fracazmix(ig) = (fraca(ig,lmix(ig)+1)-fraca(ig,lmix(ig)))/ &
703      (zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))*zmix(ig) + &
704      fraca(ig, lmix(ig)) - zlev(ig, lmix(ig))*(fraca(ig,lmix(ig)+1)-fraca(ig &
705      ,lmix(ig)))/(zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))
706  END DO
707
708  DO l = 2, nlay
709    DO ig = 1, ngrid
710      IF (larg_cons(ig,l)>1.) THEN
711        IF (l>lmix(ig)) THEN
712          ! test
713          IF (zmax(ig)-zmix(ig)<1.E-10) THEN
714            ! print*,'pb xxx'
715            xxx(ig, l) = (lmaxa(ig)+1.-l)/(lmaxa(ig)+1.-lmix(ig))
716          ELSE
717            xxx(ig, l) = (zmax(ig)-zlev(ig,l))/(zmax(ig)-zmix(ig))
718          END IF
719          IF (idetr==0) THEN
720            fraca(ig, l) = fracazmix(ig)
721          ELSE IF (idetr==1) THEN
722            fraca(ig, l) = fracazmix(ig)*xxx(ig, l)
723          ELSE IF (idetr==2) THEN
724            fraca(ig, l) = fracazmix(ig)*(1.-(1.-xxx(ig,l))**2)
725          ELSE
726            fraca(ig, l) = fracazmix(ig)*xxx(ig, l)**2
727          END IF
728          ! print*,ig,l,lmix(ig),lmaxa(ig),xxx(ig,l),'LLLLLLL'
729          fraca(ig, l) = max(fraca(ig,l), 0.)
730          fraca(ig, l) = min(fraca(ig,l), 0.5)
731          fracd(ig, l) = 1. - fraca(ig, l)
732          fracc(ig, l) = larg_cons(ig, l)/(r_aspect*zmax(ig))
733        END IF
734      END IF
735    END DO
736  END DO
737
738  ! print*,'fin calcul fraca'
739  ! print*,'11 OK convect8'
740  ! print*,'Ea3 ',wa_moy
741  ! ------------------------------------------------------------------
742  ! Calcul de fracd, wd
743  ! somme wa - wd = 0
744  ! ------------------------------------------------------------------
745
746
747  DO ig = 1, ngrid
748    fm(ig, 1) = 0.
749    fm(ig, nlay+1) = 0.
750  END DO
751
752  DO l = 2, nlay
753    DO ig = 1, ngrid
754      fm(ig, l) = fraca(ig, l)*wa_moy(ig, l)*rhobarz(ig, l)
755      ! CR:test
756      IF (entr(ig,l-1)<1E-10 .AND. fm(ig,l)>fm(ig,l-1) .AND. l>lmix(ig)) THEN
757        fm(ig, l) = fm(ig, l-1)
758        ! write(1,*)'ajustement fm, l',l
759      END IF
760      ! write(1,*)'ig,l,fm(ig,l)',ig,l,fm(ig,l)
761      ! RC
762    END DO
763    DO ig = 1, ngrid
764      IF (fracd(ig,l)<0.1) THEN
765        abort_message = 'fracd trop petit'
766        CALL abort_physic(modname, abort_message, 1)
767
768      ELSE
769        ! vitesse descendante "diagnostique"
770        wd(ig, l) = fm(ig, l)/(fracd(ig,l)*rhobarz(ig,l))
771      END IF
772    END DO
773  END DO
774
775  DO l = 1, nlay
776    DO ig = 1, ngrid
777      ! masse(ig,l)=rho(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
778      masse(ig, l) = (pplev(ig,l)-pplev(ig,l+1))/rg
779    END DO
780  END DO
781
782  ! print*,'12 OK convect8'
783  ! print*,'WA4 ',wa_moy
784  ! c------------------------------------------------------------------
785  ! calcul du transport vertical
786  ! ------------------------------------------------------------------
787
788  GO TO 4444
789  ! print*,'XXXXXXXXXXXXXXX ptimestep= ',ptimestep
790  DO l = 2, nlay - 1
791    DO ig = 1, ngrid
792      IF (fm(ig,l+1)*ptimestep>masse(ig,l) .AND. fm(ig,l+1)*ptimestep>masse( &
793          ig,l+1)) THEN
794        ! print*,'WARN!!! FM>M ig=',ig,' l=',l,'  FM='
795        ! s         ,fm(ig,l+1)*ptimestep
796        ! s         ,'   M=',masse(ig,l),masse(ig,l+1)
797      END IF
798    END DO
799  END DO
800
801  DO l = 1, nlay
802    DO ig = 1, ngrid
803      IF (entr(ig,l)*ptimestep>masse(ig,l)) THEN
804        ! print*,'WARN!!! E>M ig=',ig,' l=',l,'  E=='
805        ! s         ,entr(ig,l)*ptimestep
806        ! s         ,'   M=',masse(ig,l)
807      END IF
808    END DO
809  END DO
810
811  DO l = 1, nlay
812    DO ig = 1, ngrid
813      IF (.NOT. fm(ig,l)>=0. .OR. .NOT. fm(ig,l)<=10.) THEN
814        ! print*,'WARN!!! fm exagere ig=',ig,'   l=',l
815        ! s         ,'   FM=',fm(ig,l)
816      END IF
817      IF (.NOT. masse(ig,l)>=1.E-10 .OR. .NOT. masse(ig,l)<=1.E4) THEN
818        ! print*,'WARN!!! masse exagere ig=',ig,'   l=',l
819        ! s         ,'   M=',masse(ig,l)
820        ! print*,'rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)',
821        ! s                 rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)
822        ! print*,'zlev(ig,l+1),zlev(ig,l)'
823        ! s                ,zlev(ig,l+1),zlev(ig,l)
824        ! print*,'pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)'
825        ! s                ,pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)
826      END IF
827      IF (.NOT. entr(ig,l)>=0. .OR. .NOT. entr(ig,l)<=10.) THEN
828        ! print*,'WARN!!! entr exagere ig=',ig,'   l=',l
829        ! s         ,'   E=',entr(ig,l)
830      END IF
831    END DO
832  END DO
833
8344444 CONTINUE
835
836  ! CR:redefinition du entr
837  DO l = 1, nlay
838    DO ig = 1, ngrid
839      detr(ig, l) = fm(ig, l) + entr(ig, l) - fm(ig, l+1)
840      IF (detr(ig,l)<0.) THEN
841        ! entr(ig,l)=entr(ig,l)-detr(ig,l)
842        fm(ig, l+1) = fm(ig, l) + entr(ig, l)
843        detr(ig, l) = 0.
844        ! print*,'WARNING !!! detrainement negatif ',ig,l
845      END IF
846    END DO
847  END DO
848  ! RC
849  IF (w2di==1) THEN
850    fm0 = fm0 + ptimestep*(fm-fm0)/tho
851    entr0 = entr0 + ptimestep*(entr-entr0)/tho
852  ELSE
853    fm0 = fm
854    entr0 = entr
855  END IF
856
857  IF (1==1) THEN
858    CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zh, zdhadj, &
859      zha)
860    CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zo, pdoadj, &
861      zoa)
862  ELSE
863    CALL dqthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, zh, &
864      zdhadj, zha)
865    CALL dqthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, zo, &
866      pdoadj, zoa)
867  END IF
868
869  IF (1==0) THEN
870    CALL dvthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, zmax, &
871      zu, zv, pduadj, pdvadj, zua, zva)
872  ELSE
873    CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zu, pduadj, &
874      zua)
875    CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zv, pdvadj, &
876      zva)
877  END IF
878
879  DO l = 1, nlay
880    DO ig = 1, ngrid
881      zf = 0.5*(fracc(ig,l)+fracc(ig,l+1))
882      zf2 = zf/(1.-zf)
883      thetath2(ig, l) = zf2*(zha(ig,l)-zh(ig,l))**2
884      wth2(ig, l) = zf2*(0.5*(wa_moy(ig,l)+wa_moy(ig,l+1)))**2
885    END DO
886  END DO
887
888
889
890  ! print*,'13 OK convect8'
891  ! print*,'WA5 ',wa_moy
892  DO l = 1, nlay
893    DO ig = 1, ngrid
894      pdtadj(ig, l) = zdhadj(ig, l)*zpspsk(ig, l)
895    END DO
896  END DO
897
898
899  ! do l=1,nlay
900  ! do ig=1,ngrid
901  ! if(abs(pdtadj(ig,l))*86400..gt.500.) then
902  ! print*,'WARN!!! ig=',ig,'  l=',l
903  ! s         ,'   pdtadj=',pdtadj(ig,l)
904  ! endif
905  ! if(abs(pdoadj(ig,l))*86400..gt.1.) then
906  ! print*,'WARN!!! ig=',ig,'  l=',l
907  ! s         ,'   pdoadj=',pdoadj(ig,l)
908  ! endif
909  ! enddo
910  ! enddo
911
912  ! print*,'14 OK convect8'
913  ! ------------------------------------------------------------------
914  ! Calculs pour les sorties
915  ! ------------------------------------------------------------------
916
917  IF (sorties) THEN
918    DO l = 1, nlay
919      DO ig = 1, ngrid
920        zla(ig, l) = (1.-fracd(ig,l))*zmax(ig)
921        zld(ig, l) = fracd(ig, l)*zmax(ig)
922        IF (1.-fracd(ig,l)>1.E-10) zwa(ig, l) = wd(ig, l)*fracd(ig, l)/ &
923          (1.-fracd(ig,l))
924      END DO
925    END DO
926
927    ! deja fait
928    ! do l=1,nlay
929    ! do ig=1,ngrid
930    ! detr(ig,l)=fm(ig,l)+entr(ig,l)-fm(ig,l+1)
931    ! if (detr(ig,l).lt.0.) then
932    ! entr(ig,l)=entr(ig,l)-detr(ig,l)
933    ! detr(ig,l)=0.
934    ! print*,'WARNING !!! detrainement negatif ',ig,l
935    ! endif
936    ! enddo
937    ! enddo
938
939    ! print*,'15 OK convect8'
940
941    isplit = isplit + 1
942
943
944    ! #define und
945    GO TO 123
946#ifdef und
947    CALL writeg1d(1, nlay, wd, 'wd      ', 'wd      ')
948    CALL writeg1d(1, nlay, zwa, 'wa      ', 'wa      ')
949    CALL writeg1d(1, nlay, fracd, 'fracd      ', 'fracd      ')
950    CALL writeg1d(1, nlay, fraca, 'fraca      ', 'fraca      ')
951    CALL writeg1d(1, nlay, wa_moy, 'wam         ', 'wam         ')
952    CALL writeg1d(1, nlay, zla, 'la      ', 'la      ')
953    CALL writeg1d(1, nlay, zld, 'ld      ', 'ld      ')
954    CALL writeg1d(1, nlay, pt, 'pt      ', 'pt      ')
955    CALL writeg1d(1, nlay, zh, 'zh      ', 'zh      ')
956    CALL writeg1d(1, nlay, zha, 'zha      ', 'zha      ')
957    CALL writeg1d(1, nlay, zu, 'zu      ', 'zu      ')
958    CALL writeg1d(1, nlay, zv, 'zv      ', 'zv      ')
959    CALL writeg1d(1, nlay, zo, 'zo      ', 'zo      ')
960    CALL writeg1d(1, nlay, wh, 'wh      ', 'wh      ')
961    CALL writeg1d(1, nlay, wu, 'wu      ', 'wu      ')
962    CALL writeg1d(1, nlay, wv, 'wv      ', 'wv      ')
963    CALL writeg1d(1, nlay, wo, 'w15uo     ', 'wXo     ')
964    CALL writeg1d(1, nlay, zdhadj, 'zdhadj      ', 'zdhadj      ')
965    CALL writeg1d(1, nlay, pduadj, 'pduadj      ', 'pduadj      ')
966    CALL writeg1d(1, nlay, pdvadj, 'pdvadj      ', 'pdvadj      ')
967    CALL writeg1d(1, nlay, pdoadj, 'pdoadj      ', 'pdoadj      ')
968    CALL writeg1d(1, nlay, entr, 'entr        ', 'entr        ')
969    CALL writeg1d(1, nlay, detr, 'detr        ', 'detr        ')
970    CALL writeg1d(1, nlay, fm, 'fm          ', 'fm          ')
971
972    CALL writeg1d(1, nlay, pdtadj, 'pdtadj    ', 'pdtadj    ')
973    CALL writeg1d(1, nlay, pplay, 'pplay     ', 'pplay     ')
974    CALL writeg1d(1, nlay, pplev, 'pplev     ', 'pplev     ')
975
976    ! recalcul des flux en diagnostique...
977    ! print*,'PAS DE TEMPS ',ptimestep
978    CALL dt2f(pplev, pplay, pt, pdtadj, wh)
979    CALL writeg1d(1, nlay, wh, 'wh2     ', 'wh2     ')
980#endif
981123 CONTINUE
982
983  END IF
984
985  ! if(wa_moy(1,4).gt.1.e-10) stop
986
987  ! print*,'19 OK convect8'
988  RETURN
989END SUBROUTINE calcul_sec
990
991SUBROUTINE fermeture_seche(ngrid, nlay, pplay, pplev, pphi, zlev, rhobarz, &
992    f0, zpspsk, alim_star, zh, zo, lentr, lmin, nu_min, nu_max, r_aspect, &
993    zmax, wmax)
994
995  USE dimphy
996  IMPLICIT NONE
997
998  include "YOMCST.h"
999
1000  INTEGER ngrid, nlay
1001  REAL pplay(ngrid, nlay), pplev(ngrid, nlay+1)
1002  REAL pphi(ngrid, nlay)
1003  REAL zlev(klon, klev+1)
1004  REAL alim_star(klon, klev)
1005  REAL f0(klon)
1006  INTEGER lentr(klon)
1007  INTEGER lmin(klon)
1008  REAL zmax(klon)
1009  REAL wmax(klon)
1010  REAL nu_min
1011  REAL nu_max
1012  REAL r_aspect
1013  REAL rhobarz(klon, klev+1)
1014  REAL zh(klon, klev)
1015  REAL zo(klon, klev)
1016  REAL zpspsk(klon, klev)
1017
1018  INTEGER ig, l
1019
1020  REAL f_star(klon, klev+1)
1021  REAL detr_star(klon, klev)
1022  REAL entr_star(klon, klev)
1023  REAL zw2(klon, klev+1)
1024  REAL linter(klon)
1025  INTEGER lmix(klon)
1026  INTEGER lmax(klon)
1027  REAL zlevinter(klon)
1028  REAL wa_moy(klon, klev+1)
1029  REAL wmaxa(klon)
1030  REAL ztv(klon, klev)
1031  REAL ztva(klon, klev)
1032  REAL nu(klon, klev)
1033  ! real zmax0_sec(klon)
1034  ! save zmax0_sec
1035  REAL, SAVE, ALLOCATABLE :: zmax0_sec(:)
1036  !$OMP THREADPRIVATE(zmax0_sec)
1037  LOGICAL, SAVE :: first = .TRUE.
1038  !$OMP THREADPRIVATE(first)
1039
1040  IF (first) THEN
1041    ALLOCATE (zmax0_sec(klon))
1042    first = .FALSE.
1043  END IF
1044
1045  DO l = 1, nlay
1046    DO ig = 1, ngrid
1047      ztv(ig, l) = zh(ig, l)/zpspsk(ig, l)
1048      ztv(ig, l) = ztv(ig, l)*(1.+retv*zo(ig,l))
1049    END DO
1050  END DO
1051  DO l = 1, nlay - 2
1052    DO ig = 1, ngrid
1053      IF (ztv(ig,l)>ztv(ig,l+1) .AND. alim_star(ig,l)>1.E-10 .AND. &
1054          zw2(ig,l)<1E-10) THEN
1055        f_star(ig, l+1) = alim_star(ig, l)
1056        ! test:calcul de dteta
1057        zw2(ig, l+1) = 2.*rg*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig, l+1)* &
1058          (zlev(ig,l+1)-zlev(ig,l))*0.4*pphi(ig, l)/(pphi(ig,l+1)-pphi(ig,l))
1059      ELSE IF ((zw2(ig,l)>=1E-10) .AND. (f_star(ig,l)+alim_star(ig, &
1060          l))>1.E-10) THEN
1061        ! estimation du detrainement a partir de la geometrie du pas
1062        ! precedent
1063        ! tests sur la definition du detr
1064        nu(ig, l) = (nu_min+nu_max)/2.*(1.-(nu_max-nu_min)/(nu_max+nu_min)* &
1065          tanh((((ztva(ig,l-1)-ztv(ig,l))/ztv(ig,l))/0.0005)))
1066
1067        detr_star(ig, l) = rhobarz(ig, l)*sqrt(zw2(ig,l))/ &
1068          (r_aspect*zmax0_sec(ig))* & ! s
1069                                      ! /(r_aspect*zmax0(ig))*
1070          (sqrt(nu(ig,l)*zlev(ig,l+1)/sqrt(zw2(ig,l)))-sqrt(nu(ig,l)*zlev(ig, &
1071          l)/sqrt(zw2(ig,l))))
1072        detr_star(ig, l) = detr_star(ig, l)/f0(ig)
1073        IF ((detr_star(ig,l))>f_star(ig,l)) THEN
1074          detr_star(ig, l) = f_star(ig, l)
1075        END IF
1076        entr_star(ig, l) = 0.9*detr_star(ig, l)
1077        IF ((l<lentr(ig))) THEN
1078          entr_star(ig, l) = 0.
1079          ! detr_star(ig,l)=0.
1080        END IF
1081        ! print*,'ok detr_star'
1082        ! prise en compte du detrainement dans le calcul du flux
1083        f_star(ig, l+1) = f_star(ig, l) + alim_star(ig, l) + &
1084          entr_star(ig, l) - detr_star(ig, l)
1085        ! test sur le signe de f_star
1086        IF ((f_star(ig,l+1)+detr_star(ig,l))>1.E-10) THEN
1087          ! AM on melange Tl et qt du thermique
1088          ztva(ig, l) = (f_star(ig,l)*ztva(ig,l-1)+(entr_star(ig, &
1089            l)+alim_star(ig,l))*ztv(ig,l))/(f_star(ig,l+1)+detr_star(ig,l))
1090          zw2(ig, l+1) = zw2(ig, l)*(f_star(ig,l)/(f_star(ig, &
1091            l+1)+detr_star(ig,l)))**2 + 2.*rg*(ztva(ig,l)-ztv(ig,l))/ztv(ig, &
1092            l)*(zlev(ig,l+1)-zlev(ig,l))
1093        END IF
1094      END IF
1095
1096      IF (zw2(ig,l+1)<0.) THEN
1097        linter(ig) = (l*(zw2(ig,l+1)-zw2(ig,l))-zw2(ig,l))/(zw2(ig,l+1)-zw2( &
1098          ig,l))
1099        zw2(ig, l+1) = 0.
1100        ! print*,'linter=',linter(ig)
1101      ELSE
1102        wa_moy(ig, l+1) = sqrt(zw2(ig,l+1))
1103      END IF
1104      IF (wa_moy(ig,l+1)>wmaxa(ig)) THEN
1105        ! lmix est le niveau de la couche ou w (wa_moy) est maximum
1106        lmix(ig) = l + 1
1107        wmaxa(ig) = wa_moy(ig, l+1)
1108      END IF
1109    END DO
1110  END DO
1111  ! print*,'fin calcul zw2'
1112
1113  ! Calcul de la couche correspondant a la hauteur du thermique
1114  DO ig = 1, ngrid
1115    lmax(ig) = lentr(ig)
1116  END DO
1117  DO ig = 1, ngrid
1118    DO l = nlay, lentr(ig) + 1, -1
1119      IF (zw2(ig,l)<=1.E-10) THEN
1120        lmax(ig) = l - 1
1121      END IF
1122    END DO
1123  END DO
1124  ! pas de thermique si couche 1 stable
1125  DO ig = 1, ngrid
1126    IF (lmin(ig)>1) THEN
1127      lmax(ig) = 1
1128      lmin(ig) = 1
1129      lentr(ig) = 1
1130    END IF
1131  END DO
1132
1133  ! Determination de zw2 max
1134  DO ig = 1, ngrid
1135    wmax(ig) = 0.
1136  END DO
1137
1138  DO l = 1, nlay
1139    DO ig = 1, ngrid
1140      IF (l<=lmax(ig)) THEN
1141        IF (zw2(ig,l)<0.) THEN
1142          ! print*,'pb2 zw2<0'
1143        END IF
1144        zw2(ig, l) = sqrt(zw2(ig,l))
1145        wmax(ig) = max(wmax(ig), zw2(ig,l))
1146      ELSE
1147        zw2(ig, l) = 0.
1148      END IF
1149    END DO
1150  END DO
1151
1152  ! Longueur caracteristique correspondant a la hauteur des thermiques.
1153  DO ig = 1, ngrid
1154    zmax(ig) = 0.
1155    zlevinter(ig) = zlev(ig, 1)
1156  END DO
1157  DO ig = 1, ngrid
1158    ! calcul de zlevinter
1159    zlevinter(ig) = (zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*linter(ig) + &
1160      zlev(ig, lmax(ig)) - lmax(ig)*(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))
1161    ! pour le cas ou on prend tjs lmin=1
1162    ! zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,lmin(ig)))
1163    zmax(ig) = max(zmax(ig), zlevinter(ig)-zlev(ig,1))
1164    zmax0_sec(ig) = zmax(ig)
1165  END DO
1166
1167  RETURN
1168END SUBROUTINE fermeture_seche
Note: See TracBrowser for help on using the repository browser.