source: LMDZ6/trunk/libf/phylmdiso/thermcell_main.F90 @ 3927

Last change on this file since 3927 was 3927, checked in by Laurent Fairhead, 3 years ago

Initial import of the physics wih isotopes from Camille Risi
CR

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