source: LMDZ6/trunk/libf/phylmdiso/concvl.F90 @ 5278

Last change on this file since 5278 was 5274, checked in by abarral, 25 hours ago

Replace yomcst.h by existing module

  • Property svn:keywords set to Id
File size: 29.4 KB
Line 
1SUBROUTINE concvl(iflag_clos, &
2                  dtime, paprs, pplay, k_upper_cv, &
3                  t, q, t_wake, q_wake, s_wake, u, v, tra, ntra, &
4                  Ale, Alp, sig1, w01, &
5                  d_t, d_q, d_qcomp, d_u, d_v, d_tra, &
6                  rain, snow, kbas, ktop, sigd, &
7                  cbmf, plcl, plfc, wbeff, convoccur, &
8                  upwd, dnwd, dnwdbis, &
9                  Ma, mip, Vprecip, &
10                  cape, cin, tvp, Tconv, iflag, &
11                  pbase, bbase, dtvpdt1, dtvpdq1, dplcldt, dplcldr, &
12                  qcondc, wd, pmflxr, pmflxs, &
13!RomP >>>
14!!     .             da,phi,mp,dd_t,dd_q,lalim_conv,wght_th)
15                  da, phi, mp, phii, d1a, dam, sij, qta, clw, elij, &! RomP
16                  dd_t, dd_q, lalim_conv, wght_th, &                 ! RomP
17                  evap, ep, epmlmMm, eplaMm, &                       ! RomP
18                  wdtrainA, wdtrainS, wdtrainM, wght, qtc, sigt, detrain, &
19                  tau_cld_cv, coefw_cld_cv, &                           ! RomP+RL, AJ
20!RomP <<<
21                  epmax_diag & ! epmax_cape
22#ifdef ISO
23     &       ,xt,xt_wake,d_xt,xtrain,xtsnow &
24     &       ,xtVprecip,xtVprecipi   &
25     &       ,xtclw,dd_xt,xtevap,xtwdtrainA &
26#ifdef DIAGISO
27     &       , qlp,xtlp,qvp,xtvp & ! juste diagnostique 
28     &       , fq_detrainement,fq_ddft,fq_fluxmasse,fq_evapprecip &
29     &       , fxt_detrainement,fxt_ddft,fxt_fluxmasse,fxt_evapprecip &
30     &       , f_detrainement,q_detrainement,xt_detrainement &
31#endif         
32#endif
33     &        ) ! **************************************************************
34! *
35! CONCVL                                                      *
36! *
37! *
38! written by   : Sandrine Bony-Lena, 17/05/2003, 11.16.04    *
39! modified by :                                               *
40! **************************************************************
41
42
43  USE dimphy
44  USE infotrac_phy, ONLY: nbtr
45#ifdef ISO
46  USE infotrac_phy, ONLY: ntraciso=>ntiso
47  USE isotopes_mod, ONLY: iso_eau, bidouille_anti_divergence, ridicule, &
48        iso_eau,iso_HDO
49#endif
50#ifdef ISOVERIF
51      USE isotopes_verif_mod, ONLY: errmax,errmaxrel, &
52        iso_verif_egalite_choix,iso_verif_aberrant,iso_verif_egalite, &
53        iso_verif_noNaN,iso_verif_aberrant_encadre
54#endif
55#ifdef ISOTRAC
56      USE isotrac_routines_mod, only: iso_verif_traceur_jbid_vect
57#ifdef ISOVERIF
58      USE isotopes_verif_mod, ONLY: iso_verif_traceur_vect, &
59&       iso_verif_trac_masse_vect, iso_verif_traceur,  &
60&       iso_verif_traceur_justmass
61#endif
62#endif
63  USE phys_local_var_mod, ONLY: omega
64  USE print_control_mod, ONLY: prt_level, lunout
65  USE yomcst_mod_h, ONLY: RPI, RCLUM, RHPLA, RKBOL, RNAVO                   &
66          , RDAY, REA, REPSM, RSIYEA, RSIDAY, ROMEGA                  &
67          , R_ecc, R_peri, R_incl                                      &
68          , RA, RG, R1SA                                         &
69          , RSIGMA                                                     &
70          , R, RMD, RMV, RD, RV, RCPD                    &
71          , RMO3, RMCO2, RMC, RMCH4, RMN2O, RMCFC11, RMCFC12        &
72          , RCPV, RCVD, RCVV, RKAPPA, RETV, eps_w                    &
73          , RCW, RCS                                                 &
74          , RLVTT, RLSTT, RLMLT, RTT, RATM                           &
75          , RESTT, RALPW, RBETW, RGAMW, RALPS, RBETS, RGAMS            &
76          , RALPD, RBETD, RGAMD
77  IMPLICIT NONE
78! ======================================================================
79! Auteur(s): S. Bony-Lena (LMD/CNRS) date: ???
80! Objet: schema de convection de Emanuel (1991) interface
81! ======================================================================
82! Arguments:
83! dtime--input-R-pas d'integration (s)
84! s-------input-R-la vAleur "s" pour chaque couche
85! sigs----input-R-la vAleur "sigma" de chaque couche
86! sig-----input-R-la vAleur de "sigma" pour chaque niveau
87! psolpa--input-R-la pression au sol (en Pa)
88! pskapa--input-R-exponentiel kappa de psolpa
89! h-------input-R-enthAlpie potentielle (Cp*T/P**kappa)
90! q-------input-R-vapeur d'eau (en kg/kg)
91
92! work*: input et output: deux variables de travail,
93! on peut les mettre a 0 au debut
94! ALE--------input-R-energie disponible pour soulevement
95! ALP--------input-R-puissance disponible pour soulevement
96
97! d_h--------output-R-increment de l'enthAlpie potentielle (h)
98! d_q--------output-R-increment de la vapeur d'eau
99! rain-------output-R-la pluie (mm/s)
100! snow-------output-R-la neige (mm/s)
101! upwd-------output-R-saturated updraft mass flux (kg/m**2/s)
102! dnwd-------output-R-saturated downdraft mass flux (kg/m**2/s)
103! dnwd0------output-R-unsaturated downdraft mass flux (kg/m**2/s)
104! Ma---------output-R-adiabatic ascent mass flux (kg/m2/s)
105! mip--------output-R-mass flux shed by adiabatic ascent (kg/m2/s)
106! Vprecip----output-R-vertical profile of total precipitation (kg/m2/s)
107! Tconv------output-R-environment temperature seen by convective scheme (K)
108! Cape-------output-R-CAPE (J/kg)
109! Cin -------output-R-CIN  (J/kg)
110! Tvp--------output-R-Temperature virtuelle d'une parcelle soulevee
111! adiabatiquement a partir du niveau 1 (K)
112! deltapb----output-R-distance entre LCL et base de la colonne (<0 ; Pa)
113! Ice_flag---input-L-TRUE->prise en compte de la thermodynamique de la glace
114! dd_t-------output-R-increment de la temperature du aux descentes precipitantes
115! dd_q-------output-R-increment de la vapeur d'eau du aux desc precip
116! lalim_conv-
117! wght_th----
118! evap-------output-R
119! ep---------output-R
120! epmlmMm----output-R
121! eplaMm-----output-R
122! wdtrainA---output-R
123! wdtrainS---output-R
124! wdtrainM---output-R
125! wght-------output-R
126! ======================================================================
127
128
129  include "clesphys.h"
130
131  INTEGER, INTENT(IN)                           :: iflag_clos
132  REAL, INTENT(IN)                              :: dtime
133  REAL, DIMENSION(klon,klev),   INTENT(IN)      :: pplay
134  REAL, DIMENSION(klon,klev+1), INTENT(IN)      :: paprs
135  INTEGER,                      INTENT(IN)      :: k_upper_cv
136  REAL, DIMENSION(klon,klev),   INTENT(IN)      :: t, q, u, v
137  REAL, DIMENSION(klon,klev),   INTENT(IN)      :: t_wake, q_wake
138  REAL, DIMENSION(klon),        INTENT(IN)      :: s_wake
139  REAL, DIMENSION(klon,klev, nbtr),INTENT(IN)   :: tra
140  INTEGER,                      INTENT(IN)      :: ntra
141  REAL, DIMENSION(klon),        INTENT(IN)      :: Ale, Alp
142!CR:test: on passe lentr et alim_star des thermiques
143  INTEGER, DIMENSION(klon),     INTENT(IN)      :: lalim_conv
144  REAL, DIMENSION(klon,klev),   INTENT(IN)      :: wght_th
145#ifdef ISO
146  REAL, DIMENSION(ntraciso,klon,klev),   INTENT(IN)    ::  xt
147  REAL, DIMENSION(ntraciso,klon,klev),   INTENT(IN)    ::  xt_wake
148#endif
149
150  REAL, DIMENSION(klon,klev),   INTENT(INOUT)   :: sig1, w01
151
152  REAL, DIMENSION(klon,klev),   INTENT(OUT)     :: d_t, d_q, d_qcomp, d_u, d_v
153  REAL, DIMENSION(klon,klev, nbtr),INTENT(OUT)  :: d_tra
154  REAL, DIMENSION(klon),        INTENT(OUT)     :: rain, snow
155
156  INTEGER, DIMENSION(klon),     INTENT(OUT)     :: kbas, ktop
157  REAL, DIMENSION(klon),        INTENT(OUT)     :: sigd
158  REAL, DIMENSION(klon),        INTENT(OUT)     :: cbmf, plcl, plfc, wbeff
159  REAL, DIMENSION(klon),        INTENT(OUT)     :: convoccur
160  REAL, DIMENSION(klon,klev),   INTENT(OUT)     :: upwd, dnwd, dnwdbis
161
162!!       REAL Ma(klon,klev), mip(klon,klev),Vprecip(klon,klev)                    !jyg
163  REAL, DIMENSION(klon,klev),   INTENT(OUT)     :: Ma, mip
164  REAL, DIMENSION(klon,klev+1), INTENT(OUT)     :: Vprecip                        !jyg
165  REAL, DIMENSION(klon),        INTENT(OUT)     :: cape, cin
166  REAL, DIMENSION(klon,klev),   INTENT(OUT)     :: tvp
167  REAL, DIMENSION(klon,klev),   INTENT(OUT)     :: Tconv
168  INTEGER, DIMENSION(klon),     INTENT(OUT)     :: iflag
169  REAL, DIMENSION(klon),        INTENT(OUT)     :: pbase, bbase
170  REAL, DIMENSION(klon,klev),   INTENT(OUT)     :: dtvpdt1, dtvpdq1
171  REAL, DIMENSION(klon),        INTENT(OUT)     :: dplcldt, dplcldr
172  REAL, DIMENSION(klon,klev),   INTENT(OUT)     :: qcondc
173  REAL, DIMENSION(klon),        INTENT(OUT)     :: wd
174  REAL, DIMENSION(klon,klev+1), INTENT(OUT)     :: pmflxr, pmflxs
175
176#ifdef ISO
177  REAL, DIMENSION(ntraciso,klon,klev),   INTENT(OUT)    ::  d_xt
178  REAL, DIMENSION(ntraciso,klon),   INTENT(OUT)    ::  xtrain
179  REAL, DIMENSION(ntraciso,klon),   INTENT(OUT)    ::  xtsnow
180  REAL, DIMENSION(ntraciso,klon,klev+1),   INTENT(OUT)    ::  xtVprecip
181#endif
182
183
184  REAL, DIMENSION(klon,klev),   INTENT(OUT)     :: da, mp
185  REAL, DIMENSION(klon,klev,klev),INTENT(OUT)   :: phi
186! RomP >>>
187  REAL, DIMENSION(klon,klev,klev),INTENT(OUT)   :: phii
188  REAL, DIMENSION(klon,klev),   INTENT(OUT)     :: d1a, dam
189  REAL, DIMENSION(klon,klev,klev),INTENT(OUT)   :: sij, elij
190  REAL, DIMENSION(klon,klev),   INTENT(OUT)     :: qta
191  REAL, DIMENSION(klon,klev),   INTENT(OUT)     :: clw
192  REAL, DIMENSION(klon,klev),   INTENT(OUT)     :: dd_t, dd_q
193  REAL, DIMENSION(klon,klev),   INTENT(OUT)     :: evap, ep
194  REAL, DIMENSION(klon,klev),   INTENT(OUT)     :: eplaMm
195  REAL, DIMENSION(klon,klev,klev), INTENT(OUT)  :: epmlmMm
196  REAL, DIMENSION(klon,klev),   INTENT(OUT)     :: wdtrainA, wdtrainS, wdtrainM
197! RomP <<<
198  REAL, DIMENSION(klon,klev),   INTENT(OUT)     :: wght                       !RL
199  REAL, DIMENSION(klon,klev),   INTENT(OUT)     :: qtc
200  REAL, DIMENSION(klon,klev),   INTENT(OUT)     :: sigt, detrain
201  REAL,                         INTENT(OUT)     :: tau_cld_cv, coefw_cld_cv
202  REAL, DIMENSION(klon),        INTENT(OUT)     :: epmax_diag                ! epmax_cape
203
204#ifdef ISO
205  REAL, DIMENSION(ntraciso,klon,klev),   INTENT(OUT)     :: xtevap
206  REAL, DIMENSION(ntraciso,klon,klev),   INTENT(OUT)     :: xtwdtrainA
207  REAL, DIMENSION(ntraciso,klon,klev),   INTENT(OUT)     :: xtclw
208  REAL, DIMENSION(ntraciso,klon,klev),   INTENT(OUT)     :: dd_xt
209       ! juste diagnostique
210#ifdef DIAGISO
211  REAL, DIMENSION(klon,klev),   INTENT(OUT)              :: qlp
212  REAL, DIMENSION(ntraciso,klon,klev),   INTENT(OUT)     :: xtlp
213  REAL, DIMENSION(klon,klev),   INTENT(OUT)              :: qvp
214  REAL, DIMENSION(ntraciso,klon,klev),   INTENT(OUT)     :: xtvp
215  REAL, DIMENSION(klon,klev),   INTENT(OUT)              :: fq_detrainement
216  REAL, DIMENSION(klon,klev),   INTENT(OUT)              :: fq_fluxmasse
217  REAL, DIMENSION(klon,klev),   INTENT(OUT)              :: fq_ddft
218  REAL, DIMENSION(klon,klev),   INTENT(OUT)              :: fq_evapprecip
219  REAL, DIMENSION(ntraciso,klon,klev),   INTENT(OUT)     :: fxt_detrainement
220  REAL, DIMENSION(ntraciso,klon,klev),   INTENT(OUT)     :: fxt_fluxmasse
221  REAL, DIMENSION(ntraciso,klon,klev),   INTENT(OUT)     :: fxt_ddft
222  REAL, DIMENSION(ntraciso,klon,klev),   INTENT(OUT)     :: fxt_evapprecip
223  REAL, DIMENSION(klon,klev),   INTENT(OUT)              :: f_detrainement
224  REAL, DIMENSION(klon,klev),   INTENT(OUT)              :: q_detrainement
225  REAL, DIMENSION(ntraciso,klon,klev),   INTENT(OUT)     :: xt_detrainement
226#endif         
227#endif
228
229
230!
231!  Local
232!  ----
233  REAL, DIMENSION(klon,klev)                    :: em_p
234  REAL, DIMENSION(klon,klev+1)                  :: em_ph
235  REAL                                          :: em_sig1feed ! sigma at lower bound of feeding layer
236  REAL                                          :: em_sig2feed ! sigma at upper bound of feeding layer
237  REAL, DIMENSION(klev)                         :: em_wght ! weight density determining the feeding mixture
238  REAL, DIMENSION(klon,klev+1)                  :: Vprecipi                       !jyg
239!on enleve le save
240! SAVE em_sig1feed,em_sig2feed,em_wght
241
242  REAL, DIMENSION(klon)                         :: rflag
243  REAL, DIMENSION(klon)                         :: plim1, plim2
244  REAL, DIMENSION(klon)                         :: ptop2
245  REAL, DIMENSION(klon,klev)                    :: asupmax
246  REAL, DIMENSION(klon)                         :: supmax0, asupmaxmin
247  REAL                                          :: zx_t, zdelta, zx_qs, zcor
248!
249!   INTEGER iflag_mix
250!   SAVE iflag_mix
251  INTEGER                                       :: noff, minorig
252  INTEGER                                       :: i,j, k, itra
253  REAL, DIMENSION(klon,klev)                    :: qs, qs_wake
254!LF          SAVE cbmf
255!IM/JYG      REAL, SAVE, ALLOCATABLE :: cbmf(:)
256!!!$OMP THREADPRIVATE(cbmf)!
257  REAL, DIMENSION(klon)                         :: cbmflast
258#ifdef ISO
259REAL, DIMENSION(ntraciso,klon,klev+1)                  :: xtVprecipi
260  INTEGER                                       :: ixt
261#endif
262
263
264! Variables supplementaires liees au bilan d'energie
265! Real paire(klon)
266!LF      Real ql(klon,klev)
267! Save paire
268!LF      Save ql
269!LF      Real t1(klon,klev),q1(klon,klev)
270!LF      Save t1,q1
271! Data paire /1./
272  REAL, SAVE, ALLOCATABLE :: ql(:, :), q1(:, :), t1(:, :)
273!$OMP THREADPRIVATE(ql, q1, t1)
274        ! pas besoin d'isos ici
275
276! Variables liees au bilan d'energie et d'enthAlpi
277  REAL ztsol(klon)
278  REAL        h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot, &
279              h_qs_tot, qw_tot, ql_tot, qs_tot, ec_tot
280  SAVE        h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot, &
281              h_qs_tot, qw_tot, ql_tot, qs_tot, ec_tot
282!$OMP THREADPRIVATE(h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot)
283!$OMP THREADPRIVATE(h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot)
284  REAL        d_h_vcol, d_h_dair, d_qt, d_qw, d_ql, d_qs, d_ec
285  REAL        d_h_vcol_phy
286  REAL        fs_bound, fq_bound
287  SAVE        d_h_vcol_phy
288!$OMP THREADPRIVATE(d_h_vcol_phy)
289  REAL        zero_v(klon)
290  CHARACTER *15 ztit
291  INTEGER     ip_ebil ! PRINT level for energy conserv. diag.
292  SAVE        ip_ebil
293  DATA        ip_ebil/2/
294!$OMP THREADPRIVATE(ip_ebil)
295  INTEGER     if_ebil ! level for energy conserv. dignostics
296  SAVE        if_ebil
297  DATA        if_ebil/2/
298!$OMP THREADPRIVATE(if_ebil)
299!+jld ec_conser
300  REAL d_t_ec(klon, klev) ! tendance du a la conersion Ec -> E thermique
301  REAL zrcpd
302!-jld ec_conser
303!LF
304  INTEGER nloc
305  LOGICAL, SAVE            :: first = .TRUE.
306!$OMP THREADPRIVATE(first)
307  INTEGER, SAVE            :: itap, igout
308!$OMP THREADPRIVATE(itap, igout)
309  include "YOMCST2.h"
310  include "YOETHF.h"
311  include "FCTTRE.h"
312!jyg<
313  include "conema3.h"
314!>jyg
315
316  IF (first) THEN
317! Allocate some variables LF 04/2008
318
319!IM/JYG allocate(cbmf(klon))
320    ALLOCATE (ql(klon,klev))
321    ALLOCATE (t1(klon,klev))
322    ALLOCATE (q1(klon,klev))
323!
324    convoccur(:) = 0.
325!
326    itap = 0
327    igout = klon/2 + 1/klon
328  END IF
329! Incrementer le compteur de la physique
330  itap = itap + 1
331
332! Copy T into Tconv
333  DO k = 1, klev
334    DO i = 1, klon
335      Tconv(i, k) = t(i, k)
336    END DO
337  END DO
338
339  IF (if_ebil>=1) THEN
340    DO i = 1, klon
341      ztsol(i) = t(i, 1)
342      zero_v(i) = 0.
343      DO k = 1, klev
344        ql(i, k) = 0.
345      END DO
346    END DO
347  END IF
348
349! ym
350  snow(:) = 0
351#ifdef ISO
352      xtsnow(:,:)=0
353#endif
354
355  IF (first) THEN
356    first = .FALSE.
357
358! ===========================================================================
359! READ IN PARAMETERS FOR THE CLOSURE AND THE MIXING DISTRIBUTION
360! ===========================================================================
361
362    IF (iflag_con==3) THEN
363!      CALL cv3_inicp()
364      CALL cv3_inip()
365    END IF
366
367! ===========================================================================
368! READ IN PARAMETERS FOR CONVECTIVE INHIBITION BY TROPOS. DRYNESS
369! ===========================================================================
370
371! c$$$         open (56,file='supcrit.data')
372! c$$$         read (56,*) Supcrit1, Supcrit2
373! c$$$         close (56)
374
375    IF (prt_level>=10) WRITE (lunout, *) 'supcrit1, supcrit2', supcrit1, supcrit2
376
377! ===========================================================================
378! Initialisation pour les bilans d'eau et d'energie
379! ===========================================================================
380    IF (if_ebil>=1) d_h_vcol_phy = 0.
381
382    DO i = 1, klon
383      cbmf(i) = 0.
384!!          plcl(i) = 0.
385      sigd(i) = 0.
386    END DO
387  END IF !(first)
388
389! Initialisation a chaque pas de temps
390  plfc(:) = 0.
391  wbeff(:) = 100.
392  plcl(:) = 0.
393
394  DO k = 1, klev + 1
395    DO i = 1, klon
396      em_ph(i, k) = paprs(i, k)/100.0
397      pmflxr(i, k) = 0.
398      pmflxs(i, k) = 0.
399    END DO
400  END DO
401
402  DO k = 1, klev
403    DO i = 1, klon
404      em_p(i, k) = pplay(i, k)/100.0
405    END DO
406  END DO
407
408
409! Feeding layer
410
411  em_sig1feed = 1.
412!jyg<
413!  em_sig2feed = 0.97
414  em_sig2feed = cvl_sig2feed
415!>jyg
416! em_sig2feed = 0.8
417! Relative Weight densities
418  DO k = 1, klev
419    em_wght(k) = 1.
420  END DO
421!CRtest: couche alim des tehrmiques ponderee par a*
422! DO i = 1, klon
423! do k=1,lalim_conv(i)
424! em_wght(k)=wght_th(i,k)
425! print*,'em_wght=',em_wght(k),wght_th(i,k)
426! end do
427! END DO
428
429  IF (iflag_con==4) THEN
430    DO k = 1, klev
431      DO i = 1, klon
432        zx_t = t(i, k)
433        zdelta = max(0., sign(1.,rtt-zx_t))
434        zx_qs = min(0.5, r2es*foeew(zx_t,zdelta)/em_p(i,k)/100.0)
435        zcor = 1./(1.-retv*zx_qs)
436        qs(i, k) = zx_qs*zcor
437      END DO
438      DO i = 1, klon
439        zx_t = t_wake(i, k)
440        zdelta = max(0., sign(1.,rtt-zx_t))
441        zx_qs = min(0.5, r2es*foeew(zx_t,zdelta)/em_p(i,k)/100.0)
442        zcor = 1./(1.-retv*zx_qs)
443        qs_wake(i, k) = zx_qs*zcor
444      END DO
445    END DO
446  ELSE ! iflag_con=3 (modif de puristes qui fait la diffce pour la convergence numerique)
447    DO k = 1, klev
448      DO i = 1, klon
449        zx_t = t(i, k)
450        zdelta = max(0., sign(1.,rtt-zx_t))
451        zx_qs = r2es*foeew(zx_t, zdelta)/em_p(i, k)/100.0
452        zx_qs = min(0.5, zx_qs)
453        zcor = 1./(1.-retv*zx_qs)
454        zx_qs = zx_qs*zcor
455        qs(i, k) = zx_qs
456      END DO
457      DO i = 1, klon
458        zx_t = t_wake(i, k)
459        zdelta = max(0., sign(1.,rtt-zx_t))
460        zx_qs = r2es*foeew(zx_t, zdelta)/em_p(i, k)/100.0
461        zx_qs = min(0.5, zx_qs)
462        zcor = 1./(1.-retv*zx_qs)
463        zx_qs = zx_qs*zcor
464        qs_wake(i, k) = zx_qs
465      END DO
466    END DO
467  END IF ! iflag_con
468
469! ------------------------------------------------------------------
470
471! Main driver for convection:
472!                   iflag_con=3 -> nvlle version de KE (JYG)
473!                   iflag_con = 30  -> equivAlent to convect3
474!                   iflag_con = 4  -> equivAlent to convect1/2
475
476
477  IF (iflag_con==30) THEN
478
479 
480#ifdef ISO         
481#ifdef ISOVERIF
482       do k = 1, klev
483        do i = 1, klon               
484         do ixt=1,ntraciso         
485             call iso_verif_noNaN(xt(ixt,i,k),'concvl 394')
486         enddo
487        enddo !do i = 1, klon
488       enddo !do k = 1, klev       
489       if (iso_eau.gt.0) then
490       do k = 1, klev
491        do i = 1, klon   
492          call iso_verif_egalite_choix(xt(iso_eau,i,k),q(i,k), &
493     &            'concvl 174',errmax,errmaxrel)
494        enddo !do i = 1, klon
495       enddo !do k = 1, klev   
496       endif !if (iso_eau.gt.0) then
497       if (iso_HDO.gt.0) then
498       do k = 1, klev
499        do i = 1, klon         
500         if (q(i,k).gt.ridicule) then
501          call iso_verif_aberrant(xt(iso_hdo,i,k)/q(i,k),'concvl 175')
502         endif ! if (q(i,k).gt.ridicule) then
503        enddo
504       enddo   
505       endif !if (iso_eau.gt.0) then
506#ifdef ISOTRAC
507        do k = 1, klev
508        do i = 1, klon   
509           call iso_verif_traceur(xt(1,i,k),'concvl 218')
510        enddo
511        enddo
512#endif       
513       write(*,*) 'concvl 170: avant appel cv_driver'
514#endif
515        ! ISOVERIF ! end verif       
516#endif
517
518! print *, '-> cv_driver'      !jyg
519    CALL cv_driver(klon, klev, klevp1, ntra, iflag_con, &
520                   t, q, qs, u, v, tra, &
521                   em_p, em_ph, iflag, &
522                   d_t, d_q, d_u, d_v, d_tra, rain, &
523                   Vprecip, cbmf, sig1, w01, & !jyg
524                   kbas, ktop, &
525                   dtime, Ma, upwd, dnwd, dnwdbis, qcondc, wd, cape, &
526                   da, phi, mp, phii, d1a, dam, sij, clw, elij, &       !RomP
527                   evap, ep, epmlmMm, eplaMm, &                         !RomP
528                   wdtrainA, wdtrainM, &                                !RomP
529                   epmax_diag & ! epmax_cape
530#ifdef ISO
531     &             , xt,d_xt &
532     &             , xtrain,xtVprecip &
533     &             , xtevap,xtclw,xtwdtrainA &
534#ifdef DIAGISO
535     &          , qlp,xtlp,qvp,xtvp & ! juste diagnostique 
536     &          , fq_detrainement,fq_ddft,fq_fluxmasse,fq_evapprecip &
537     &          , fxt_detrainement,fxt_ddft,fxt_fluxmasse,fxt_evapprecip &
538     &          , f_detrainement,q_detrainement,xt_detrainement &
539#endif     
540#endif
541     &              )
542!           print *, 'cv_driver ->'      !jyg
543
544#ifdef ISO
545      ! verif
546#ifdef ISOVERIF
547       write(*,*) 'concvl 463: après appel cv_driver'
548       do k = 1, klev
549        do i = 1, klon
550        if (iso_eau.gt.0) then
551            call iso_verif_egalite(xt(iso_eau,i,k),q(i,k),'concvl 203')
552            call iso_verif_egalite(xt_wake(iso_eau,i,k),q_wake(i,k),'concvl 204')
553            call iso_verif_egalite(d_xt(iso_eau,i,k),d_q(i,k), &
554     &            'concvl 452')
555         endif !if (iso_eau.gt.0) then
556#ifdef DIAGISO
557         do ixt=1,ntraciso
558            call iso_verif_noNaN(xt(ixt,i,k),'concvl 460')
559            call iso_verif_noNaN(xtlp(ixt,i,k),'concvl 295')
560            call iso_verif_noNaN(xtvp(ixt,i,k),'concvl 260')
561          enddo
562#endif                 
563        enddo
564       enddo       
565#ifdef ISOTRAC
566           call iso_verif_traceur_vect(xt,klon,klev,'concvl 218')
567           call iso_verif_trac_masse_vect(d_xt,klon,klev, &
568     &           'concvl 464',errmax,errmaxrel)
569#endif           
570#endif
571       ! end verif       
572#endif
573
574    DO i = 1, klon
575      cbmf(i) = Ma(i, kbas(i))
576    END DO
577
578!RL
579    wght(:, :) = 0.
580    DO i = 1, klon
581      wght(i, 1) = 1.
582    END DO
583!RL
584
585  ELSE
586
587!LF   necessary for gathered fields
588    nloc = klon
589#ifdef ISOVERIF
590        write(*,*) 'concvl 581: juste avant appel de cva_driver'
591#endif
592    CALL cva_driver(klon, klev, klev+1, ntra, nloc, k_upper_cv, &
593                    iflag_con, iflag_mix, iflag_ice_thermo, &
594                    iflag_clos, ok_conserv_q, dtime, cvl_comp_threshold, &
595                    t, q, qs, t_wake, q_wake, qs_wake, s_wake, u, v, tra, &
596                    em_p, em_ph, &
597                    Ale, Alp, omega, &
598                    em_sig1feed, em_sig2feed, em_wght, &
599                    iflag, d_t, d_q, d_qcomp, d_u, d_v, d_tra, rain, kbas, ktop, &
600                    cbmf, plcl, plfc, wbeff, sig1, w01, ptop2, sigd, &
601                    Ma, mip, Vprecip, Vprecipi, upwd, dnwd, dnwdbis, qcondc, wd, &
602                    cape, cin, tvp, &
603                    dd_t, dd_q, plim1, plim2, asupmax, supmax0, &
604                    asupmaxmin, lalim_conv, &
605!AC!+!RomP+jyg
606!!                   da,phi,mp,phii,d1a,dam,sij,clw,elij, &               ! RomP
607!!                   evap,ep,epmlmMm,eplaMm,                              ! RomP
608                    da, phi, mp, phii, d1a, dam, sij, wght, &           ! RomP+RL
609                    qta, clw, elij, evap, ep, epmlmMm, eplaMm, &        ! RomP+RL
610                    wdtrainA, wdtrainS, wdtrainM, qtc, sigt, detrain, &
611                    tau_cld_cv, coefw_cld_cv, &                         ! RomP,AJ
612!AC!+!RomP+jyg
613                    epmax_diag & ! epmax_cape
614#ifdef ISO
615     &             ,xt,xt_wake,d_xt, xtrain &
616     &             ,xtvprecip,xtvprecipi &
617     &             ,xtclw,dd_xt,xtevap,xtwdtrainA &
618#ifdef DIAGISO     
619     &          , qlp,xtlp,qvp,xtvp &
620     &          , fq_detrainement,fq_ddft,fq_fluxmasse,fq_evapprecip &
621     &          , fxt_detrainement,fxt_ddft,fxt_fluxmasse,fxt_evapprecip &
622     &          , f_detrainement,q_detrainement,xt_detrainement &
623#endif     
624#endif
625     &          )     
626  END IF
627! ------------------------------------------------------------------
628  IF (prt_level>=10) WRITE (lunout, *) ' cva_driver -> cbmf,plcl,plfc,wbeff, d_t, d_q ', &
629                                         cbmf(1), plcl(1), plfc(1), wbeff(1), d_t(1,1), d_q(1,1)
630
631  DO i = 1, klon
632    rain(i) = rain(i)/86400.
633    rflag(i) = iflag(i)
634#ifdef ISO
635       do ixt = 1, ntraciso
636        xtrain(ixt,i) = xtrain(ixt,i)/86400.
637       enddo
638#endif
639  END DO
640
641  DO k = 1, klev
642    DO i = 1, klon
643      d_t(i, k) = dtime*d_t(i, k)
644      d_q(i, k) = dtime*d_q(i, k)
645      d_u(i, k) = dtime*d_u(i, k)
646      d_v(i, k) = dtime*d_v(i, k)
647#ifdef ISO
648           do ixt = 1, ntraciso
649            d_xt(ixt,i,k) = dtime*d_xt(ixt,i,k)
650           enddo
651#endif
652    END DO
653  END DO
654
655             
656#ifdef ISO
657#ifdef ISOVERIF     
658!        k=1
659!        i=  538
660        write(*,*) 'concvl 640'
661!        write(*,*) 'q(i,k),d_q(i,k)=', q(i,k),d_q(i,k)
662!        write(*,*) 'xt(iso_HDO,i,k),d_xt(iso_HDO,i,k)=', &
663!     &          xt(iso_HDO,i,k),d_xt(iso_HDO,i,k)
664  DO k = 1, klev
665    DO i = 1, klon
666           if (iso_HDO.gt.0) then
667             if (q(i,k).gt.ridicule) then
668                 call iso_verif_aberrant_encadre((xt(iso_HDO,i,k) &
669     &        +d_xt(iso_HDO,i,k))/(q(i,k)+d_q(i,k)),'concvl 250')
670             endif !if (q_seri(i,k).gt.ridicule) then
671          endif !if (iso_HDO.gt.0) then
672           if (iso_eau.gt.0) then
673             call iso_verif_egalite_choix(d_xt(iso_eau,i,k), &
674     &          d_q(i,k),'concvl 530',errmax*dtime,errmaxrel)
675          endif !if (iso_HDO.gt.0) then
676#ifdef ISOTRAC
677           call iso_verif_traceur_justmass(d_xt(1,i,k),'concvl 316')
678#endif 
679    END DO
680  END DO         
681#endif           
682#endif
683
684#ifdef ISO
685      if ((iso_eau.gt.0).and.(bidouille_anti_divergence)) then
686        do k=1,klev   
687        do i=1,klon
688            d_xt(iso_eau,i,k)=d_q(i,k)
689        enddo !do i=1,klon
690        enddo !do k=1,klev               
691#ifdef ISOTRAC 
692        call iso_verif_traceur_jbid_vect(d_xt, &
693     &            klon,klev)   
694#endif         
695      endif !if ((iso_eau.gt.0).and.(bidouille_anti_divergence)) then 
696#endif
697
698  IF (iflag_con==3) THEN
699    DO i = 1,klon
700      IF (wbeff(i) > 100. .OR. wbeff(i) == 0 .OR. iflag(i) > 3) THEN
701        wbeff(i) = 0.
702        convoccur(i) = 0. 
703      ELSE
704        convoccur(i) = 1.
705      ENDIF
706    ENDDO
707  ENDIF
708
709  IF (iflag_con==30) THEN
710    DO itra = 1, ntra
711      DO k = 1, klev
712        DO i = 1, klon
713!RL!            d_tra(i,k,itra) =dtime*d_tra(i,k,itra)
714          d_tra(i, k, itra) = 0.
715        END DO
716      END DO
717    END DO
718  END IF
719
720!!AC!
721  IF (iflag_con==3) THEN
722    DO itra = 1, ntra
723      DO k = 1, klev
724        DO i = 1, klon
725!RL!            d_tra(i,k,itra) =dtime*d_tra(i,k,itra)
726          d_tra(i, k, itra) = 0.
727        END DO
728      END DO
729    END DO
730  END IF
731!!AC!
732
733  DO k = 1, klev
734    DO i = 1, klon
735      t1(i, k) = t(i, k) + d_t(i, k)
736      q1(i, k) = q(i, k) + d_q(i, k)
737! juste diag: pas besoin d'isos ici
738    END DO
739  END DO
740
741!                                                     !jyg
742  IF (iflag_con == 30 .OR. iflag_ice_thermo ==0) THEN
743! --Separation neige/pluie (pour diagnostics)         !jyg
744    DO k = 1, klev                                    !jyg
745      DO i = 1, klon                                  !jyg
746        IF (t1(i,k)<rtt) THEN                         !jyg
747          pmflxs(i, k) = Vprecip(i, k)                !jyg
748        ELSE                                          !jyg
749          pmflxr(i, k) = Vprecip(i, k)                !jyg
750        END IF                                        !jyg
751      END DO                                          !jyg
752    END DO                                            !jyg
753  ELSE
754    DO k = 1, klev                                    !jyg
755      DO i = 1, klon                                  !jyg
756        pmflxs(i, k) = Vprecipi(i, k)                 !jyg
757        pmflxr(i, k) = Vprecip(i, k)-Vprecipi(i, k)   !jyg
758      END DO                                          !jyg
759    END DO                                            !jyg
760  ENDIF
761
762! c      IF (if_ebil.ge.2) THEN
763! c        ztit='after convect'
764! c        CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
765! c     e      , t1,q1,ql,qs,u,v,paprs,pplay
766! c     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
767! c         call diagphy(paire,ztit,ip_ebil
768! c     e      , zero_v, zero_v, zero_v, zero_v, zero_v
769! c     e      , zero_v, rain, zero_v, ztsol
770! c     e      , d_h_vcol, d_qt, d_ec
771! c     s      , fs_bound, fq_bound )
772! c      END IF
773
774
775! les traceurs ne sont pas mis dans cette version de convect4:
776  IF (iflag_con==4) THEN
777    DO itra = 1, ntra
778      DO k = 1, klev
779        DO i = 1, klon
780          d_tra(i, k, itra) = 0.
781        END DO
782      END DO
783    END DO
784  END IF
785! print*, 'concvl->: dd_t,dd_q ',dd_t(1,1),dd_q(1,1)
786
787  DO k = 1, klev
788    DO i = 1, klon
789      dtvpdt1(i, k) = 0.
790      dtvpdq1(i, k) = 0.
791    END DO
792  END DO
793  DO i = 1, klon
794    dplcldt(i) = 0.
795    dplcldr(i) = 0.
796  END DO
797
798  IF (prt_level>=20) THEN
799    DO k = 1, klev
800! print*,'physiq apres_add_con i k it d_u d_v d_t d_q qdl0',igout, &
801!         k,itap,d_u_con(igout,k) ,d_v_con(igout,k), d_t_con(igout,k), &
802!         d_q_con(igout,k),dql0(igout,k)
803! print*,'phys apres_add_con itap Ma cin ALE ALP wak t q undi t q', &
804!         itap,Ma(igout,k),cin(igout),ALE(igout), ALP(igout), &
805!         t_wake(igout,k),q_wake(igout,k),t_undi(igout,k),q_undi(igout,k)
806! print*,'phy apres_add_con itap CON rain snow EMA wk1 wk2 Vpp mip', &
807!         itap,rain_con(igout),snow_con(igout),ema_work1(igout,k), &
808!         ema_work2(igout,k),Vprecip(igout,k), mip(igout,k)
809! print*,'phy apres_add_con itap upwd dnwd dnwd0 cape tvp Tconv ', &
810!         itap,upwd(igout,k),dnwd(igout,k),dnwd0(igout,k),cape(igout), &
811!         tvp(igout,k),Tconv(igout,k)
812! print*,'phy apres_add_con itap dtvpdt dtvdq dplcl dplcldr qcondc', &
813!         itap,dtvpdt1(igout,k),dtvpdq1(igout,k),dplcldt(igout), &
814!         dplcldr(igout),qcondc(igout,k)
815! print*,'phy apres_add_con itap wd pmflxr Kpmflxr Kp1 Kpmflxs Kp1', &
816!         itap,wd(igout),pmflxr(igout,k),pmflxr(igout,k+1),pmflxs(igout,k), &
817!         pmflxs(igout,k+1)
818! print*,'phy apres_add_con itap da phi mp ftd fqd lalim wgth', &
819!         itap,da(igout,k),phi(igout,k,k),mp(igout,k),ftd(igout,k), &
820!         fqd(igout,k),lalim_conv(igout),wght_th(igout,k)
821    END DO
822  END IF !(prt_level.EQ.20) THEN
823
824  RETURN
825END SUBROUTINE concvl
826
Note: See TracBrowser for help on using the repository browser.