source: LMDZ6/trunk/libf/phylmd/thermcell_main.F90 @ 3897

Last change on this file since 3897 was 3451, checked in by fhourdin, 6 years ago

Saving old version of thermcell_plume in thermcell_plume_6A.F90
It includes both thermcell_plume_6A and thermcell_plume_5B corresponding
to the 5B and 6A versions used for CMIP5 and CMIP6.
The latest was previously named thermcellV1_plume.
The new thermcell_plume is a clean version (removing obsolete
options) of thermcell_plume_6A.
The 3 versions are controled by
flag_thermals_ed <= 9 thermcell_plume_6A

<= 19 thermcell_plume_5B
else thermcell_plume (default 20 for convergence with 6A)

Fredho

  • 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:keywords set to Author Date Id Revision
File size: 32.7 KB
Line 
1!
2! $Id: thermcell_main.F90 3451 2019-01-27 11:07:30Z jyg $
3!
4      SUBROUTINE thermcell_main(itap,ngrid,nlay,ptimestep  &
5     &                  ,pplay,pplev,pphi,debut  &
6     &                  ,pu,pv,pt,po  &
7     &                  ,pduadj,pdvadj,pdtadj,pdoadj  &
8     &                  ,fm0,entr0,detr0,zqta,zqla,lmax  &
9     &                  ,ratqscth,ratqsdiff,zqsatth  &
10     &                  ,Ale_bl,Alp_bl,lalim_conv,wght_th &
11     &                  ,zmax0, f0,zw2,fraca,ztv &
12     &                  ,zpspsk,ztla,zthl &
13!!! nrlmd le 10/04/2012
14     &                  ,pbl_tke,pctsrf,omega,airephy &
15     &                  ,zlcl,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
16     &                  ,n2,s2,ale_bl_stat &
17     &                  ,therm_tke_max,env_tke_max &
18     &                  ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
19     &                  ,alp_bl_conv,alp_bl_stat &
20!!! fin nrlmd le 10/04/2012
21     &                  ,ztva  )
22
23      USE dimphy
24      USE ioipsl
25      USE indice_sol_mod
26      USE print_control_mod, ONLY: lunout,prt_level
27      IMPLICIT NONE
28
29!=======================================================================
30!   Auteurs: Frederic Hourdin, Catherine Rio, Anne Mathieu
31!   Version du 09.02.07
32!   Calcul du transport vertical dans la couche limite en presence
33!   de "thermiques" explicitement representes avec processus nuageux
34!
35!   Reecriture a partir d'un listing papier a Habas, le 14/02/00
36!
37!   le thermique est suppose homogene et dissipe par melange avec
38!   son environnement. la longueur l_mix controle l'efficacite du
39!   melange
40!
41!   Le calcul du transport des differentes especes se fait en prenant
42!   en compte:
43!     1. un flux de masse montant
44!     2. un flux de masse descendant
45!     3. un entrainement
46!     4. un detrainement
47!
48! Modif 2013/01/04 (FH hourdin@lmd.jussieu.fr)
49!    Introduction of an implicit computation of vertical advection in
50!    the environment of thermal plumes in thermcell_dq
51!    impl =     0 : explicit, 1 : implicit, -1 : old version
52!    controled by iflag_thermals =
53!       15, 16 run with impl=-1 : numerical convergence with NPv3
54!       17, 18 run with impl=1  : more stable
55!    15 and 17 correspond to the activation of the stratocumulus "bidouille"
56!
57!=======================================================================
58
59
60!-----------------------------------------------------------------------
61!   declarations:
62!   -------------
63
64#include "YOMCST.h"
65#include "YOETHF.h"
66#include "FCTTRE.h"
67#include "thermcell.h"
68
69!   arguments:
70!   ----------
71
72!IM 140508
73      INTEGER itap
74
75      INTEGER ngrid,nlay
76      real ptimestep
77      REAL pt(ngrid,nlay),pdtadj(ngrid,nlay)
78      REAL pu(ngrid,nlay),pduadj(ngrid,nlay)
79      REAL pv(ngrid,nlay),pdvadj(ngrid,nlay)
80      REAL po(ngrid,nlay),pdoadj(ngrid,nlay)
81      REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
82      real pphi(ngrid,nlay)
83      LOGICAL debut
84
85!   local:
86!   ------
87
88      integer icount
89
90      integer, save :: dvdq=1,dqimpl=-1
91!$OMP THREADPRIVATE(dvdq,dqimpl)
92      data icount/0/
93      save icount
94!$OMP THREADPRIVATE(icount)
95
96      integer,save :: igout=1
97!$OMP THREADPRIVATE(igout)
98      integer,save :: lunout1=6
99!$OMP THREADPRIVATE(lunout1)
100      integer,save :: lev_out=10
101!$OMP THREADPRIVATE(lev_out)
102
103      REAL susqr2pi, Reuler
104
105      INTEGER ig,k,l,ll,ierr
106      real zsortie1d(klon)
107      INTEGER lmax(klon),lmin(klon),lalim(klon)
108      INTEGER lmix(klon)
109      INTEGER lmix_bis(klon)
110      real linter(klon)
111      real zmix(klon)
112      real zmax(klon),zw2(klon,klev+1),ztva(klon,klev),zw_est(klon,klev+1),ztva_est(klon,klev)
113!      real fraca(klon,klev)
114
115      real zmax_sec(klon)
116!on garde le zmax du pas de temps precedent
117      real zmax0(klon)
118!FH/IM     save zmax0
119
120      real lambda
121
122      real zlev(klon,klev+1),zlay(klon,klev)
123      real deltaz(klon,klev)
124      REAL zh(klon,klev)
125      real zthl(klon,klev),zdthladj(klon,klev)
126      REAL ztv(klon,klev)
127      real zu(klon,klev),zv(klon,klev),zo(klon,klev)
128      real zl(klon,klev)
129      real zsortie(klon,klev)
130      real zva(klon,klev)
131      real zua(klon,klev)
132      real zoa(klon,klev)
133
134      real zta(klon,klev)
135      real zha(klon,klev)
136      real fraca(klon,klev+1)
137      real zf,zf2
138      real thetath2(klon,klev),wth2(klon,klev),wth3(klon,klev)
139      real q2(klon,klev)
140! FH probleme de dimensionnement avec l'allocation dynamique
141!     common/comtherm/thetath2,wth2
142      real wq(klon,klev)
143      real wthl(klon,klev)
144      real wthv(klon,klev)
145   
146      real ratqscth(klon,klev)
147      real var
148      real vardiff
149      real ratqsdiff(klon,klev)
150
151      logical sorties
152      real rho(klon,klev),rhobarz(klon,klev),masse(klon,klev)
153      real zpspsk(klon,klev)
154
155      real wmax(klon)
156      real wmax_tmp(klon)
157      real wmax_sec(klon)
158      real fm0(klon,klev+1),entr0(klon,klev),detr0(klon,klev)
159      real fm(klon,klev+1),entr(klon,klev),detr(klon,klev)
160
161      real ztla(klon,klev),zqla(klon,klev),zqta(klon,klev)
162!niveau de condensation
163      integer nivcon(klon)
164      real zcon(klon)
165      REAL CHI
166      real zcon2(klon)
167      real pcon(klon)
168      real zqsat(klon,klev)
169      real zqsatth(klon,klev)
170
171      real f_star(klon,klev+1),entr_star(klon,klev)
172      real detr_star(klon,klev)
173      real alim_star_tot(klon)
174      real alim_star(klon,klev)
175      real alim_star_clos(klon,klev)
176      real f(klon), f0(klon)
177!FH/IM     save f0
178      real zlevinter(klon)
179       real seuil
180      real csc(klon,klev)
181
182!!! nrlmd le 10/04/2012
183
184!------Entrées
185      real pbl_tke(klon,klev+1,nbsrf)
186      real pctsrf(klon,nbsrf)
187      real omega(klon,klev)
188      real airephy(klon)
189!------Sorties
190      real zlcl(klon),fraca0(klon),w0(klon),w_conv(klon)
191      real therm_tke_max0(klon),env_tke_max0(klon)
192      real n2(klon),s2(klon)
193      real ale_bl_stat(klon)
194      real therm_tke_max(klon,klev),env_tke_max(klon,klev)
195      real alp_bl_det(klon),alp_bl_fluct_m(klon),alp_bl_fluct_tke(klon),alp_bl_conv(klon),alp_bl_stat(klon)
196!------Local
197      integer nsrf
198      real rhobarz0(klon)                    ! Densité au LCL
199      logical ok_lcl(klon)                   ! Existence du LCL des thermiques
200      integer klcl(klon)                     ! Niveau du LCL
201      real interp(klon)                      ! Coef d'interpolation pour le LCL
202!--Triggering
203      real Su                                ! Surface unité: celle d'un updraft élémentaire
204      parameter(Su=4e4)
205      real hcoef                             ! Coefficient directeur pour le calcul de s2
206      parameter(hcoef=1)
207      real hmincoef                          ! Coefficient directeur pour l'ordonnée à l'origine pour le calcul de s2
208      parameter(hmincoef=0.3)
209      real eps1                              ! Fraction de surface occupée par la population 1 : eps1=n1*s1/(fraca0*Sd)
210      parameter(eps1=0.3)
211      real hmin(ngrid)                       ! Ordonnée à l'origine pour le calcul de s2
212      real zmax_moy(ngrid)                   ! Hauteur moyenne des thermiques : zmax_moy = zlcl + 0.33 (zmax-zlcl)
213      real zmax_moy_coef
214      parameter(zmax_moy_coef=0.33)
215      real depth(klon)                       ! Epaisseur moyenne du cumulus
216      real w_max(klon)                       ! Vitesse max statistique
217      real s_max(klon)
218!--Closure
219      real pbl_tke_max(klon,klev)            ! Profil de TKE moyenne
220      real pbl_tke_max0(klon)                ! TKE moyenne au LCL
221      real w_ls(klon,klev)                   ! Vitesse verticale grande échelle (m/s)
222      real coef_m                            ! On considère un rendement pour alp_bl_fluct_m
223      parameter(coef_m=1.)
224      real coef_tke                          ! On considère un rendement pour alp_bl_fluct_tke
225      parameter(coef_tke=1.)
226
227!!! fin nrlmd le 10/04/2012
228
229!
230      !nouvelles variables pour la convection
231      real Ale_bl(klon)
232      real Alp_bl(klon)
233      real alp_int(klon),dp_int(klon),zdp
234      real ale_int(klon)
235      integer n_int(klon)
236      real fm_tot(klon)
237      real wght_th(klon,klev)
238      integer lalim_conv(klon)
239!v1d     logical therm
240!v1d     save therm
241
242      character*2 str2
243      character*10 str10
244
245      character (len=20) :: modname='thermcell_main'
246      character (len=80) :: abort_message
247
248      EXTERNAL SCOPY
249!
250
251!-----------------------------------------------------------------------
252!   initialisation:
253!   ---------------
254!
255
256   seuil=0.25
257
258   if (debut) then
259      if (iflag_thermals==15.or.iflag_thermals==16) then
260         dvdq=0
261         dqimpl=-1
262      else
263         dvdq=1
264         dqimpl=1
265      endif
266
267      fm0=0.
268      entr0=0.
269      detr0=0.
270   endif
271   fm=0. ; entr=0. ; detr=0.
272   icount=icount+1
273
274!IM 090508 beg
275!print*,'====================================================================='
276!print*,'====================================================================='
277!print*,' PAS ',icount,' PAS ',icount,' PAS ',icount,' PAS ',icount
278!print*,'====================================================================='
279!print*,'====================================================================='
280!IM 090508 end
281
282      if (prt_level.ge.1) print*,'thermcell_main V4'
283
284       sorties=.true.
285      IF(ngrid.NE.klon) THEN
286         PRINT*
287         PRINT*,'STOP dans convadj'
288         PRINT*,'ngrid    =',ngrid
289         PRINT*,'klon  =',klon
290      ENDIF
291!
292!     write(lunout,*)'WARNING thermcell_main f0=max(f0,1.e-2)'
293     do ig=1,klon
294         f0(ig)=max(f0(ig),1.e-2)
295         zmax0(ig)=max(zmax0(ig),40.)
296!IMmarche pas ?!       if (f0(ig)<1.e-2) f0(ig)=1.e-2
297     enddo
298
299      if (prt_level.ge.20) then
300       do ig=1,ngrid
301          print*,'th_main ig f0',ig,f0(ig)
302       enddo
303      endif
304!-----------------------------------------------------------------------
305! Calcul de T,q,ql a partir de Tl et qT dans l environnement
306!   --------------------------------------------------------------------
307!
308      CALL thermcell_env(ngrid,nlay,po,pt,pu,pv,pplay,  &
309     &           pplev,zo,zh,zl,ztv,zthl,zu,zv,zpspsk,zqsat,lev_out)
310       
311      if (prt_level.ge.1) print*,'thermcell_main apres thermcell_env'
312
313!------------------------------------------------------------------------
314!                       --------------------
315!
316!
317!                       + + + + + + + + + + +
318!
319!
320!  wa, fraca, wd, fracd --------------------   zlev(2), rhobarz
321!  wh,wt,wo ...
322!
323!                       + + + + + + + + + + +  zh,zu,zv,zo,rho
324!
325!
326!                       --------------------   zlev(1)
327!                       \\\\\\\\\\\\\\\\\\\\
328!
329!
330
331!-----------------------------------------------------------------------
332!   Calcul des altitudes des couches
333!-----------------------------------------------------------------------
334
335      do l=2,nlay
336         zlev(:,l)=0.5*(pphi(:,l)+pphi(:,l-1))/RG
337      enddo
338         zlev(:,1)=0.
339         zlev(:,nlay+1)=(2.*pphi(:,klev)-pphi(:,klev-1))/RG
340      do l=1,nlay
341         zlay(:,l)=pphi(:,l)/RG
342      enddo
343!calcul de l epaisseur des couches
344      do l=1,nlay
345         deltaz(:,l)=zlev(:,l+1)-zlev(:,l)
346      enddo
347
348!     print*,'2 OK convect8'
349!-----------------------------------------------------------------------
350!   Calcul des densites
351!-----------------------------------------------------------------------
352
353     rho(:,:)=pplay(:,:)/(zpspsk(:,:)*RD*ztv(:,:))
354
355     if (prt_level.ge.10)write(lunout,*)                                &
356    &    'WARNING thermcell_main rhobarz(:,1)=rho(:,1)'
357      rhobarz(:,1)=rho(:,1)
358
359      do l=2,nlay
360         rhobarz(:,l)=0.5*(rho(:,l)+rho(:,l-1))
361      enddo
362
363!calcul de la masse
364      do l=1,nlay
365         masse(:,l)=(pplev(:,l)-pplev(:,l+1))/RG
366      enddo
367
368      if (prt_level.ge.1) print*,'thermcell_main apres initialisation'
369
370!------------------------------------------------------------------
371!
372!             /|\
373!    --------  |  F_k+1 -------   
374!                              ----> D_k
375!             /|\              <---- E_k , A_k
376!    --------  |  F_k ---------
377!                              ----> D_k-1
378!                              <---- E_k-1 , A_k-1
379!
380!
381!
382!
383!
384!    ---------------------------
385!
386!    ----- F_lmax+1=0 ----------         \
387!            lmax     (zmax)              |
388!    ---------------------------          |
389!                                         |
390!    ---------------------------          |
391!                                         |
392!    ---------------------------          |
393!                                         |
394!    ---------------------------          |
395!                                         |
396!    ---------------------------          |
397!                                         |  E
398!    ---------------------------          |  D
399!                                         |
400!    ---------------------------          |
401!                                         |
402!    ---------------------------  \       |
403!            lalim                 |      |
404!    ---------------------------   |      |
405!                                  |      |
406!    ---------------------------   |      |
407!                                  | A    |
408!    ---------------------------   |      |
409!                                  |      |
410!    ---------------------------   |      |
411!    lmin  (=1 pour le moment)     |      |
412!    ----- F_lmin=0 ------------  /      /
413!
414!    ---------------------------
415!    //////////////////////////
416!
417!
418!=============================================================================
419!  Calculs initiaux ne faisant pas intervenir les changements de phase
420!=============================================================================
421
422!------------------------------------------------------------------
423!  1. alim_star est le profil vertical de l'alimentation a la base du
424!     panache thermique, calcule a partir de la flotabilite de l'air sec
425!  2. lmin et lalim sont les indices inferieurs et superieurs de alim_star
426!------------------------------------------------------------------
427!
428      entr_star=0. ; detr_star=0. ; alim_star=0. ; alim_star_tot=0.
429      lmin=1
430
431!-----------------------------------------------------------------------------
432!  3. wmax_sec et zmax_sec sont les vitesses et altitudes maximum d'un
433!     panache sec conservatif (e=d=0) alimente selon alim_star
434!     Il s'agit d'un calcul de type CAPE
435!     zmax_sec est utilise pour determiner la geometrie du thermique.
436!------------------------------------------------------------------------------
437!---------------------------------------------------------------------------------
438!calcul du melange et des variables dans le thermique
439!--------------------------------------------------------------------------------
440!
441      if (prt_level.ge.1) print*,'avant thermcell_plume ',lev_out
442
443!=====================================================================
444! Old version of thermcell_plume in thermcell_plume_6A.F90
445! It includes both thermcell_plume_6A and thermcell_plume_5B corresponding
446! to the 5B and 6A versions used for CMIP5 and CMIP6.
447! The latest was previously named thermcellV1_plume.
448! The new thermcell_plume is a clean version (removing obsolete
449! options) of thermcell_plume_6A.
450! The 3 versions are controled by
451! flag_thermals_ed <= 9 thermcell_plume_6A
452!                  <= 19 thermcell_plume_5B
453!                  else thermcell_plume (default 20 for convergence with 6A)
454! Fredho
455!=====================================================================
456
457      if (iflag_thermals_ed<=9) then
458!         print*,'THERM NOUVELLE/NOUVELLE Arnaud'
459         CALL thermcell_plume_6A(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,&
460     &    zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot,  &
461     &    lalim,f0,detr_star,entr_star,f_star,csc,ztva,  &
462     &    ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
463     &    ,lev_out,lunout1,igout)
464
465      elseif (iflag_thermals_ed<=19) then
466!        print*,'THERM RIO et al 2010, version d Arnaud'
467         CALL thermcell_plume_5B(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,&
468     &    zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot,  &
469     &    lalim,f0,detr_star,entr_star,f_star,csc,ztva,  &
470     &    ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
471     &    ,lev_out,lunout1,igout)
472      else
473         CALL thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,&
474     &    zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot,  &
475     &    lalim,f0,detr_star,entr_star,f_star,csc,ztva,  &
476     &    ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
477     &    ,lev_out,lunout1,igout)
478      endif
479
480      if (prt_level.ge.1) print*,'apres thermcell_plume ',lev_out
481
482      call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lalim ')
483      call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lmix  ')
484
485      if (prt_level.ge.1) print*,'thermcell_main apres thermcell_plume'
486      if (prt_level.ge.10) then
487         write(lunout1,*) 'Dans thermcell_main 2'
488         write(lunout1,*) 'lmin ',lmin(igout)
489         write(lunout1,*) 'lalim ',lalim(igout)
490         write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
491         write(lunout1,'(i6,i4,4e15.5)') (igout,l,alim_star(igout,l),entr_star(igout,l),detr_star(igout,l) &
492     &    ,f_star(igout,l+1),l=1,nint(linter(igout))+5)
493      endif
494
495!-------------------------------------------------------------------------------
496! Calcul des caracteristiques du thermique:zmax,zmix,wmax
497!-------------------------------------------------------------------------------
498!
499      CALL thermcell_height(ngrid,nlay,lalim,lmin,linter,lmix,zw2,  &
500     &           zlev,lmax,zmax,zmax0,zmix,wmax,lev_out)
501! Attention, w2 est transforme en sa racine carree dans cette routine
502! Le probleme vient du fait que linter et lmix sont souvent égaux à 1.
503      wmax_tmp=0.
504      do  l=1,nlay
505         wmax_tmp(:)=max(wmax_tmp(:),zw2(:,l))
506      enddo
507!     print*,"ZMAX ",lalim,lmin,linter,lmix,lmax,zmax,zmax0,zmix,wmax
508
509
510
511      call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lalim ')
512      call test_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmin  ')
513      call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmix  ')
514      call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmax  ')
515
516      if (prt_level.ge.1) print*,'thermcell_main apres thermcell_height'
517
518!-------------------------------------------------------------------------------
519! Fermeture,determination de f
520!-------------------------------------------------------------------------------
521!
522!
523!!      write(lunout,*)'THERM NOUVEAU XXXXX'
524      CALL thermcell_dry(ngrid,nlay,zlev,pphi,ztv,alim_star,  &
525    &                      lalim,lmin,zmax_sec,wmax_sec,lev_out)
526
527 
528call test_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry  lmin  ')
529call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry  lalim ')
530
531      if (prt_level.ge.1) print*,'thermcell_main apres thermcell_dry'
532      if (prt_level.ge.10) then
533         write(lunout1,*) 'Dans thermcell_main 1b'
534         write(lunout1,*) 'lmin ',lmin(igout)
535         write(lunout1,*) 'lalim ',lalim(igout)
536         write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
537         write(lunout1,'(i6,i4,e15.5)') (igout,l,alim_star(igout,l) &
538     &    ,l=1,lalim(igout)+4)
539      endif
540
541
542
543
544! Choix de la fonction d'alimentation utilisee pour la fermeture.
545! Apparemment sans importance
546      alim_star_clos(:,:)=alim_star(:,:)
547      alim_star_clos(:,:)=entr_star(:,:)+alim_star(:,:)
548!
549!CR Appel de la fermeture seche
550      if (iflag_thermals_closure.eq.1) then
551
552      CALL thermcell_closure(ngrid,nlay,r_aspect_thermals,ptimestep,rho,  &
553     &   zlev,lalim,alim_star_clos,f_star,zmax_sec,wmax_sec,f,lev_out)
554
555!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
556! Appel avec les zmax et wmax tenant compte de la condensation
557! Semble moins bien marcher
558     else if (iflag_thermals_closure.eq.2) then
559
560     CALL thermcell_closure(ngrid,nlay,r_aspect_thermals,ptimestep,rho,  &
561    &   zlev,lalim,alim_star,f_star,zmax,wmax,f,lev_out)
562
563     endif
564
565!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
566
567      if(prt_level.ge.1)print*,'thermcell_closure apres thermcell_closure'
568
569      if (tau_thermals>1.) then
570         lambda=exp(-ptimestep/tau_thermals)
571         f0=(1.-lambda)*f+lambda*f0
572      else
573         f0=f
574      endif
575
576! Test valable seulement en 1D mais pas genant
577      if (.not. (f0(1).ge.0.) ) then
578              abort_message = '.not. (f0(1).ge.0.)'
579              CALL abort_physic (modname,abort_message,1)
580      endif
581
582!-------------------------------------------------------------------------------
583!deduction des flux
584!-------------------------------------------------------------------------------
585
586      CALL thermcell_flux2(ngrid,nlay,ptimestep,masse, &
587     &       lalim,lmax,alim_star,  &
588     &       entr_star,detr_star,f,rhobarz,zlev,zw2,fm,entr,  &
589     &       detr,zqla,lev_out,lunout1,igout)
590!IM 060508    &       detr,zqla,zmax,lev_out,lunout,igout)
591
592      if (prt_level.ge.1) print*,'thermcell_main apres thermcell_flux'
593      call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lalim ')
594      call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lmax  ')
595
596!------------------------------------------------------------------
597!   On ne prend pas directement les profils issus des calculs precedents
598!   mais on s'autorise genereusement une relaxation vers ceci avec
599!   une constante de temps tau_thermals (typiquement 1800s).
600!------------------------------------------------------------------
601
602      if (tau_thermals>1.) then
603         lambda=exp(-ptimestep/tau_thermals)
604         fm0=(1.-lambda)*fm+lambda*fm0
605         entr0=(1.-lambda)*entr+lambda*entr0
606         detr0=(1.-lambda)*detr+lambda*detr0
607      else
608         fm0=fm
609         entr0=entr
610         detr0=detr
611      endif
612
613!c------------------------------------------------------------------
614!   calcul du transport vertical
615!------------------------------------------------------------------
616
617      call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse,  &
618     &                    zthl,zdthladj,zta,lev_out)
619      call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse,  &
620     &                   po,pdoadj,zoa,lev_out)
621
622!------------------------------------------------------------------
623! Calcul de la fraction de l'ascendance
624!------------------------------------------------------------------
625      do ig=1,klon
626         fraca(ig,1)=0.
627         fraca(ig,nlay+1)=0.
628      enddo
629      do l=2,nlay
630         do ig=1,klon
631            if (zw2(ig,l).gt.1.e-10) then
632            fraca(ig,l)=fm(ig,l)/(rhobarz(ig,l)*zw2(ig,l))
633            else
634            fraca(ig,l)=0.
635            endif
636         enddo
637      enddo
638     
639!------------------------------------------------------------------
640!  calcul du transport vertical du moment horizontal
641!------------------------------------------------------------------
642
643!IM 090508 
644      if (dvdq == 0 ) then
645
646! Calcul du transport de V tenant compte d'echange par gradient
647! de pression horizontal avec l'environnement
648
649         call thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,masse  &
650!    &    ,fraca*dvdq,zmax &
651     &    ,fraca,zmax &
652     &    ,zu,zv,pduadj,pdvadj,zua,zva,lev_out)
653
654      else
655
656! calcul purement conservatif pour le transport de V
657         call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse  &
658     &    ,zu,pduadj,zua,lev_out)
659         call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse  &
660     &    ,zv,pdvadj,zva,lev_out)
661
662      endif
663
664!     print*,'13 OK convect8'
665      do l=1,nlay
666         do ig=1,ngrid
667           pdtadj(ig,l)=zdthladj(ig,l)*zpspsk(ig,l) 
668         enddo
669      enddo
670
671      if (prt_level.ge.1) print*,'14 OK convect8'
672!------------------------------------------------------------------
673!   Calculs de diagnostiques pour les sorties
674!------------------------------------------------------------------
675!calcul de fraca pour les sorties
676     
677      if (sorties) then
678      if (prt_level.ge.1) print*,'14a OK convect8'
679! calcul du niveau de condensation
680! initialisation
681      do ig=1,ngrid
682         nivcon(ig)=0
683         zcon(ig)=0.
684      enddo
685!nouveau calcul
686      do ig=1,ngrid
687      CHI=zh(ig,1)/(1669.0-122.0*zo(ig,1)/zqsat(ig,1)-zh(ig,1))
688      pcon(ig)=pplay(ig,1)*(zo(ig,1)/zqsat(ig,1))**CHI
689      enddo
690!IM   do k=1,nlay
691      do k=1,nlay-1
692         do ig=1,ngrid
693         if ((pcon(ig).le.pplay(ig,k))  &
694     &      .and.(pcon(ig).gt.pplay(ig,k+1))) then
695            zcon2(ig)=zlay(ig,k)-(pcon(ig)-pplay(ig,k))/(RG*rho(ig,k))/100.
696         endif
697         enddo
698      enddo
699!IM
700      ierr=0
701      do ig=1,ngrid
702        if (pcon(ig).le.pplay(ig,nlay)) then
703           zcon2(ig)=zlay(ig,nlay)-(pcon(ig)-pplay(ig,nlay))/(RG*rho(ig,nlay))/100.
704           ierr=1
705        endif
706      enddo
707      if (ierr==1) then
708           abort_message = 'thermcellV0_main: les thermiques vont trop haut '
709           CALL abort_physic (modname,abort_message,1)
710      endif
711
712      if (prt_level.ge.1) print*,'14b OK convect8'
713      do k=nlay,1,-1
714         do ig=1,ngrid
715            if (zqla(ig,k).gt.1e-10) then
716               nivcon(ig)=k
717               zcon(ig)=zlev(ig,k)
718            endif
719         enddo
720      enddo
721      if (prt_level.ge.1) print*,'14c OK convect8'
722!calcul des moments
723!initialisation
724      do l=1,nlay
725         do ig=1,ngrid
726            q2(ig,l)=0.
727            wth2(ig,l)=0.
728            wth3(ig,l)=0.
729            ratqscth(ig,l)=0.
730            ratqsdiff(ig,l)=0.
731         enddo
732      enddo     
733      if (prt_level.ge.1) print*,'14d OK convect8'
734      if (prt_level.ge.10)write(lunout,*)                                &
735    &     'WARNING thermcell_main wth2=0. si zw2 > 1.e-10'
736      do l=1,nlay
737         do ig=1,ngrid
738            zf=fraca(ig,l)
739            zf2=zf/(1.-zf)
740!
741            thetath2(ig,l)=zf2*(ztla(ig,l)-zthl(ig,l))**2
742            if(zw2(ig,l).gt.1.e-10) then
743             wth2(ig,l)=zf2*(zw2(ig,l))**2
744            else
745             wth2(ig,l)=0.
746            endif
747            wth3(ig,l)=zf2*(1-2.*fraca(ig,l))/(1-fraca(ig,l))  &
748     &                *zw2(ig,l)*zw2(ig,l)*zw2(ig,l)
749            q2(ig,l)=zf2*(zqta(ig,l)*1000.-po(ig,l)*1000.)**2
750!test: on calcul q2/po=ratqsc
751            ratqscth(ig,l)=sqrt(max(q2(ig,l),1.e-6)/(po(ig,l)*1000.))
752         enddo
753      enddo
754!calcul des flux: q, thetal et thetav
755      do l=1,nlay
756         do ig=1,ngrid
757      wq(ig,l)=fraca(ig,l)*zw2(ig,l)*(zqta(ig,l)*1000.-po(ig,l)*1000.)
758      wthl(ig,l)=fraca(ig,l)*zw2(ig,l)*(ztla(ig,l)-zthl(ig,l))
759      wthv(ig,l)=fraca(ig,l)*zw2(ig,l)*(ztva(ig,l)-ztv(ig,l))
760         enddo
761      enddo
762!
763! $Id: thermcell_main.F90 3451 2019-01-27 11:07:30Z jyg $
764!
765      CALL thermcell_alp(ngrid,nlay,ptimestep  &
766     &                  ,pplay,pplev  &
767     &                  ,fm0,entr0,lmax  &
768     &                  ,Ale_bl,Alp_bl,lalim_conv,wght_th &
769     &                  ,zw2,fraca &
770!!! necessire en plus
771     &                  ,pcon,rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax &
772!!! nrlmd le 10/04/2012
773     &                  ,pbl_tke,pctsrf,omega,airephy &
774     &                  ,zlcl,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
775     &                  ,n2,s2,ale_bl_stat &
776     &                  ,therm_tke_max,env_tke_max &
777     &                  ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
778     &                  ,alp_bl_conv,alp_bl_stat &
779!!! fin nrlmd le 10/04/2012
780     &                   )
781
782
783
784!calcul du ratqscdiff
785      if (prt_level.ge.1) print*,'14e OK convect8'
786      var=0.
787      vardiff=0.
788      ratqsdiff(:,:)=0.
789
790      do l=1,klev
791         do ig=1,ngrid
792            if (l<=lalim(ig)) then
793            var=var+alim_star(ig,l)*zqta(ig,l)*1000.
794            endif
795         enddo
796      enddo
797
798      if (prt_level.ge.1) print*,'14f OK convect8'
799
800      do l=1,klev
801         do ig=1,ngrid
802            if (l<=lalim(ig)) then
803               zf=fraca(ig,l)
804               zf2=zf/(1.-zf)
805               vardiff=vardiff+alim_star(ig,l)*(zqta(ig,l)*1000.-var)**2
806            endif
807         enddo
808      enddo
809
810      if (prt_level.ge.1) print*,'14g OK convect8'
811      do l=1,nlay
812         do ig=1,ngrid
813            ratqsdiff(ig,l)=sqrt(vardiff)/(po(ig,l)*1000.)   
814!           write(11,*)'ratqsdiff=',ratqsdiff(ig,l)
815         enddo
816      enddo
817!--------------------------------------------------------------------   
818!
819!ecriture des fichiers sortie
820!     print*,'15 OK convect8 CCCCCCCCCCCCCCCCCCc'
821
822      endif
823
824      if (prt_level.ge.1) print*,'thermcell_main FIN  OK'
825
826      return
827      end
828
829!-----------------------------------------------------------------------------
830
831      subroutine test_ltherm(klon,klev,pplev,pplay,long,seuil,ztv,po,ztva,zqla,f_star,zw2,comment)
832      USE print_control_mod, ONLY: prt_level
833      IMPLICIT NONE
834
835      integer i, k, klon,klev
836      real pplev(klon,klev+1),pplay(klon,klev)
837      real ztv(klon,klev)
838      real po(klon,klev)
839      real ztva(klon,klev)
840      real zqla(klon,klev)
841      real f_star(klon,klev)
842      real zw2(klon,klev)
843      integer long(klon)
844      real seuil
845      character*21 comment
846
847      if (prt_level.ge.1) THEN
848       print*,'WARNING !!! TEST ',comment
849      endif
850      return
851
852!  test sur la hauteur des thermiques ...
853         do i=1,klon
854!IMtemp           if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
855           if (prt_level.ge.10) then
856               print*,'WARNING ',comment,' au point ',i,' K= ',long(i)
857               print*,'  K  P(MB)  THV(K)     Qenv(g/kg)THVA        QLA(g/kg)   F*        W2'
858               do k=1,klev
859                  write(6,'(i3,7f10.3)') k,pplay(i,k),ztv(i,k),1000*po(i,k),ztva(i,k),1000*zqla(i,k),f_star(i,k),zw2(i,k)
860               enddo
861           endif
862         enddo
863
864
865      return
866      end
867
868!!! nrlmd le 10/04/2012                          Transport de la TKE par le thermique moyen pour la fermeture en ALP
869!                                                         On transporte pbl_tke pour donner therm_tke
870!                                          Copie conforme de la subroutine DTKE dans physiq.F écrite par Frederic Hourdin
871      subroutine thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0,  &
872     &           rg,pplev,therm_tke_max)
873      USE print_control_mod, ONLY: prt_level
874      implicit none
875
876!=======================================================================
877!
878!   Calcul du transport verticale dans la couche limite en presence
879!   de "thermiques" explicitement representes
880!   calcul du dq/dt une fois qu'on connait les ascendances
881!
882!=======================================================================
883
884      integer ngrid,nlay,nsrf
885
886      real ptimestep
887      real masse0(ngrid,nlay),fm0(ngrid,nlay+1),pplev(ngrid,nlay+1)
888      real entr0(ngrid,nlay),rg
889      real therm_tke_max(ngrid,nlay)
890      real detr0(ngrid,nlay)
891
892
893      real masse(ngrid,nlay),fm(ngrid,nlay+1)
894      real entr(ngrid,nlay)
895      real q(ngrid,nlay)
896      integer lev_out                           ! niveau pour les print
897
898      real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
899
900      real zzm
901
902      integer ig,k
903      integer isrf
904
905
906      lev_out=0
907
908
909      if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
910
911!   calcul du detrainement
912      do k=1,nlay
913         detr0(:,k)=fm0(:,k)-fm0(:,k+1)+entr0(:,k)
914         masse0(:,k)=(pplev(:,k)-pplev(:,k+1))/RG
915      enddo
916
917
918! Decalage vertical des entrainements et detrainements.
919      masse(:,1)=0.5*masse0(:,1)
920      entr(:,1)=0.5*entr0(:,1)
921      detr(:,1)=0.5*detr0(:,1)
922      fm(:,1)=0.
923      do k=1,nlay-1
924         masse(:,k+1)=0.5*(masse0(:,k)+masse0(:,k+1))
925         entr(:,k+1)=0.5*(entr0(:,k)+entr0(:,k+1))
926         detr(:,k+1)=0.5*(detr0(:,k)+detr0(:,k+1))
927         fm(:,k+1)=fm(:,k)+entr(:,k)-detr(:,k)
928      enddo
929      fm(:,nlay+1)=0.
930
931!!! nrlmd le 16/09/2010
932!   calcul de la valeur dans les ascendances
933!       do ig=1,ngrid
934!          qa(ig,1)=q(ig,1)
935!       enddo
936!!!
937
938!do isrf=1,nsrf
939
940!   q(:,:)=therm_tke(:,:,isrf)
941   q(:,:)=therm_tke_max(:,:)
942!!! nrlmd le 16/09/2010
943      do ig=1,ngrid
944         qa(ig,1)=q(ig,1)
945      enddo
946!!!
947
948    if (1==1) then
949      do k=2,nlay
950         do ig=1,ngrid
951            if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.  &
952     &         1.e-5*masse(ig,k)) then
953         qa(ig,k)=(fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k))  &
954     &         /(fm(ig,k+1)+detr(ig,k))
955            else
956               qa(ig,k)=q(ig,k)
957            endif
958            if (qa(ig,k).lt.0.) then
959!               print*,'qa<0!!!'
960            endif
961            if (q(ig,k).lt.0.) then
962!               print*,'q<0!!!'
963            endif
964         enddo
965      enddo
966
967! Calcul du flux subsident
968
969      do k=2,nlay
970         do ig=1,ngrid
971            wqd(ig,k)=fm(ig,k)*q(ig,k)
972            if (wqd(ig,k).lt.0.) then
973!               print*,'wqd<0!!!'
974            endif
975         enddo
976      enddo
977      do ig=1,ngrid
978         wqd(ig,1)=0.
979         wqd(ig,nlay+1)=0.
980      enddo
981
982! Calcul des tendances
983      do k=1,nlay
984         do ig=1,ngrid
985            q(ig,k)=q(ig,k)+(detr(ig,k)*qa(ig,k)-entr(ig,k)*q(ig,k)  &
986     &               -wqd(ig,k)+wqd(ig,k+1))  &
987     &               *ptimestep/masse(ig,k)
988         enddo
989      enddo
990
991 endif
992
993   therm_tke_max(:,:)=q(:,:)
994
995      return
996!!! fin nrlmd le 10/04/2012
997     end
998
Note: See TracBrowser for help on using the repository browser.