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

Last change on this file since 2051 was 1992, checked in by lguez, 10 years ago

Converted to free source form files in libf/phylmd which were still in
fixed source form. The conversion was done using the polish mode of
the NAG Fortran Compiler.

In addition to converting to free source form, the processing of the
files also:

-- indented the code (including comments);

-- set Fortran keywords to uppercase, and set all other identifiers
to lower case;

-- added qualifiers to end statements (for example "end subroutine
conflx", instead of "end");

-- changed the terminating statements of all DO loops so that each
loop ends with an ENDDO statement (instead of a labeled continue).

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