LMDZ
thermcell.F90
Go to the documentation of this file.
1 
2 ! $Id: thermcell.F90 2346 2015-08-21 15:13:46Z emillour $
3 
4 SUBROUTINE 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 
834 4444 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
981 123 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
989 END SUBROUTINE calcul_sec
990 
991 SUBROUTINE 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
1168 END SUBROUTINE fermeture_seche
dimphy::klon
integer, save klon
Definition: dimphy.F90:3
dimphy::klev
integer, save klev
Definition: dimphy.F90:7
scopy
subroutine scopy(n, sx, incx, sy, incy)
Definition: cray.F:9
dvthermcell2
subroutine dvthermcell2(ngrid, nlay, ptimestep, fm, entr, masse, fraca, larga, u, v, du, dv, ua, va)
Definition: thermcell_old.F90:4395
false
!$Id itapm1 ENDIF!IM on interpole les champs sur les niveaux STD de pression!IM a chaque pas de temps de la physique c!positionnement de l argument logique a false c!pour ne pas recalculer deux fois la meme chose!c!a cet effet un appel a plevel_new a ete deplace c!a la fin de la serie d appels c!la boucle DO nlevSTD a ete internalisee c!dans d ou la creation de cette routine c c!CALL false
Definition: calcul_STDlev.h:26
calcul_sec
subroutine calcul_sec(ngrid, nlay, ptimestep, pplay, pplev, pphi, zlev, pu,pv, pt, po, zmax, wmax, zw2, lmix
Definition: thermcell.F90:7
parameter
!$Header!integer nvarmx parameter(nfmx=10, imx=200, jmx=150, lmx=200, nvarmx=1000) real xd(imx
true
!$Id itapm1 ENDIF!IM on interpole les champs sur les niveaux STD de pression!IM a chaque pas de temps de la physique c!positionnement de l argument logique a false c!pour ne pas recalculer deux fois la meme chose!c!a cet effet un appel a plevel_new a ete deplace c!a la fin de la serie d appels c!la boucle DO nlevSTD a ete internalisee c!dans d ou la creation de cette routine c c!CALL ulevSTD CALL &zphi philevSTD CALL &zx_rh rhlevSTD!DO klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon du jour ou toutes les read_climoz CALL true
Definition: calcul_STDlev.h:139
dqthermcell
subroutine dqthermcell(ngrid, nlay, ptimestep, fm, entr, masse, q, dq, qa)
Definition: thermcell_old.F90:4123
abort_physic
subroutine abort_physic(modname, message, ierr)
Definition: abort_physic.F90:3
pt
INTERFACE SUBROUTINE RRTM_ECRT_140GP pt
Definition: rrtm_ecrt_140gp.intfb.h:3
fermeture_seche
subroutine fermeture_seche(ngrid, nlay, pplay, pplev, pphi, zlev, rhobarz, f0, zpspsk, alim_star, zh, zo, lentr, lmin, nu_min, nu_max, r_aspect, zmax, wmax)
Definition: thermcell.F90:994
dimphy
Definition: dimphy.F90:1
rg
real rg
Definition: comcstphy.h:1
dqthermcell2
subroutine dqthermcell2(ngrid, nlay, ptimestep, fm, entr, masse, frac, q, dq, qa)
Definition: thermcell_old.F90:4320