source: LMDZ6/trunk/libf/phylmdiso/lmdz_wake.F90 @ 5285

Last change on this file since 5285 was 5285, checked in by abarral, 4 days ago

As discussed internally, remove generic ONLY: ... for new _mod_h modules
File size: 81.9 KB
Line 
1MODULE lmdz_wake
2
3! $Id: lmdz_wake.F90 4908 2024-04-15 17:30:55Z jyg $
4
5CONTAINS
6SUBROUTINE wake(znatsurf, p, ph, pi, dtime, &
7                te0, qe0, omgb, &
8                dtdwn, dqdwn, amdwn, amup, dta, dqa, wgen, &
9                sigd_con, Cin, &
10                deltatw, deltaqw, sigmaw, awdens, wdens, &                  ! state variables
11                dth, hw, wape, fip, gfl, &
12                dtls, dqls, ktopw, omgbdth, dp_omgb, tu, qu, &
13                dtke, dqke, omg, dp_deltomg, spread, cstar, &
14                d_deltat_gw, &
15                d_deltatw2, d_deltaqw2, d_sigmaw2, d_awdens2, d_wdens2 &     ! tendencies
16
17#ifdef ISO
18                     ,xte0,dxtdwn,dxta,deltaxtw &
19                     ,dxtls,xtu,dxtke,d_deltaxtw2 &
20#endif                 
21                ) 
22
23
24  ! **************************************************************
25  ! *
26  ! WAKE                                                        *
27  ! retour a un Pupper fixe                                *
28  ! *
29  ! written by   :  GRANDPEIX Jean-Yves   09/03/2000            *
30  ! modified by :   ROEHRIG Romain        01/29/2007            *
31  ! **************************************************************
32
33  USE ioipsl_getin_p_mod, ONLY : getin_p
34  USE dimphy
35  use mod_phys_lmdz_para
36  USE print_control_mod, ONLY: prt_level
37#ifdef ISO
38  USE infotrac_phy, ONLY : ntraciso=>ntiso
39#ifdef ISOVERIF
40  USE isotopes_verif_mod
41!, ONLY: errmax,errmaxrel
42  USE isotopes_mod, ONLY: iso_eau,iso_hdo,ridicule
43#endif
44#endif
45  USE yomcst_mod_h
46  USE cvthermo_mod_h
47  IMPLICIT NONE
48  ! ============================================================================
49
50
51  ! But : Decrire le comportement des poches froides apparaissant dans les
52  ! grands systemes convectifs, et fournir l'energie disponible pour
53  ! le declenchement de nouvelles colonnes convectives.
54
55  ! State variables :
56  ! deltatw    : temperature difference between wake and off-wake regions
57  ! deltaqw    : specific humidity difference between wake and off-wake regions
58  ! sigmaw     : fractional area covered by wakes.
59  ! wdens      : number of wakes per unit area
60
61  ! Variable de sortie :
62
63  ! wape : WAke Potential Energy
64  ! fip  : Front Incident Power (W/m2) - ALP
65  ! gfl  : Gust Front Length per unit area (m-1)
66  ! dtls : large scale temperature tendency due to wake
67  ! dqls : large scale humidity tendency due to wake
68  ! hw   : wake top hight (given by hw*deltatw(1)/2=wape)
69  ! dp_omgb : vertical gradient of large scale omega
70  ! awdens  : densite de poches actives
71  ! wdens   : densite de poches
72  ! omgbdth: flux of Delta_Theta transported by LS omega
73  ! dtKE   : differential heating (wake - unpertubed)
74  ! dqKE   : differential moistening (wake - unpertubed)
75  ! omg    : Delta_omg =vertical velocity diff. wake-undist. (Pa/s)
76  ! dp_deltomg  : vertical gradient of omg (s-1)
77  ! spread  : spreading term in d_t_wake and d_q_wake
78  ! deltatw     : updated temperature difference (T_w-T_u).
79  ! deltaqw     : updated humidity difference (q_w-q_u).
80  ! sigmaw      : updated wake fractional area.
81  ! d_deltat_gw : delta T tendency due to GW
82
83  ! Variables d'entree :
84
85  ! aire : aire de la maille
86  ! te0  : temperature dans l'environnement  (K)
87  ! qe0  : humidite dans l'environnement     (kg/kg)
88  ! omgb : vitesse verticale moyenne sur la maille (Pa/s)
89  ! dtdwn: source de chaleur due aux descentes (K/s)
90  ! dqdwn: source d'humidite due aux descentes (kg/kg/s)
91  ! dta  : source de chaleur due courants satures et detrain  (K/s)
92  ! dqa  : source d'humidite due aux courants satures et detra (kg/kg/s)
93  ! wgen : number of wakes generated per unit area and per sec (/m^2/s)
94  ! amdwn: flux de masse total des descentes, par unite de
95  !        surface de la maille (kg/m2/s)
96  ! amup : flux de masse total des ascendances, par unite de
97  !        surface de la maille (kg/m2/s)
98  ! sigd_con:
99  ! Cin  : convective inhibition
100  ! p    : pressions aux milieux des couches (Pa)
101  ! ph   : pressions aux interfaces (Pa)
102  ! pi  : (p/p_0)**kapa (adim)
103  ! dtime: increment temporel (s)
104
105  ! Variables internes :
106
107  ! rhow : masse volumique de la poche froide
108  ! rho  : environment density at P levels
109  ! rhoh : environment density at Ph levels
110  ! te   : environment temperature | may change within
111  ! qe   : environment humidity    | sub-time-stepping
112  ! the  : environment potential temperature
113  ! thu  : potential temperature in undisturbed area
114  ! tu   :  temperature  in undisturbed area
115  ! qu   : humidity in undisturbed area
116  ! dp_omgb: vertical gradient og LS omega
117  ! omgbw  : wake average vertical omega
118  ! dp_omgbw: vertical gradient of omgbw
119  ! omgbdq : flux of Delta_q transported by LS omega
120  ! dth  : potential temperature diff. wake-undist.
121  ! th1  : first pot. temp. for vertical advection (=thu)
122  ! th2  : second pot. temp. for vertical advection (=thw)
123  ! q1   : first humidity for vertical advection
124  ! q2   : second humidity for vertical advection
125  ! d_deltatw   : terme de redistribution pour deltatw
126  ! d_deltaqw   : terme de redistribution pour deltaqw
127  ! deltatw0   : deltatw initial
128  ! deltaqw0   : deltaqw initial
129  ! hw0    : wake top hight (defined as the altitude at which deltatw=0)
130  ! amflux : horizontal mass flux through wake boundary
131  ! wdens_ref: initial number of wakes per unit area (3D) or per
132  ! unit length (2D), at the beginning of each time step
133  ! Tgw    : 1 sur la période de onde de gravité
134  ! Cgw    : vitesse de propagation de onde de gravité
135  ! LL     : distance entre 2 poches
136
137  ! Arguments en entree
138  ! --------------------
139
140  INTEGER, DIMENSION (klon),        INTENT(IN)          :: znatsurf
141  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: p, pi
142  REAL, DIMENSION (klon, klev+1),   INTENT(IN)          :: ph
143  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: omgb
144  REAL,                             INTENT(IN)          :: dtime
145  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: te0, qe0
146  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: dtdwn, dqdwn
147  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: amdwn, amup
148  REAL, DIMENSION (klon, klev),     INTENT(IN)          :: dta, dqa
149  REAL, DIMENSION (klon),           INTENT(IN)          :: wgen
150  REAL, DIMENSION (klon),           INTENT(IN)          :: sigd_con
151  REAL, DIMENSION (klon),           INTENT(IN)          :: Cin
152#ifdef ISO
153  REAL, DIMENSION (ntraciso,klon, klev),     INTENT(IN)          :: xte0
154  REAL, DIMENSION (ntraciso,klon, klev),     INTENT(IN)          :: dxtdwn
155  REAL, DIMENSION (ntraciso,klon, klev),     INTENT(IN)          :: dxta
156#endif
157
158  !
159  ! Input/Output
160  ! State variables
161  REAL, DIMENSION (klon, klev),     INTENT(INOUT)       :: deltatw, deltaqw
162  REAL, DIMENSION (klon),           INTENT(INOUT)       :: sigmaw
163  REAL, DIMENSION (klon),           INTENT(INOUT)       :: awdens
164  REAL, DIMENSION (klon),           INTENT(INOUT)       :: wdens
165#ifdef ISO
166  REAL, DIMENSION (ntraciso,klon, klev),     INTENT(INOUT)       :: deltaxtw
167#endif
168
169  ! Sorties
170  ! --------
171
172  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dth
173  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: tu, qu
174  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dtls, dqls
175  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dtke, dqke
176  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: spread    !  unused (jyg)
177  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: omgbdth, omg
178  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: dp_omgb, dp_deltomg
179  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: d_deltat_gw
180  REAL, DIMENSION (klon),           INTENT(OUT)         :: hw, wape, fip, gfl, cstar
181  INTEGER, DIMENSION (klon),        INTENT(OUT)         :: ktopw
182  ! Tendencies of state variables
183  REAL, DIMENSION (klon, klev),     INTENT(OUT)         :: d_deltatw2, d_deltaqw2
184  REAL, DIMENSION (klon),           INTENT(OUT)         :: d_sigmaw2, d_awdens2, d_wdens2
185#ifdef ISO
186  REAL, DIMENSION (ntraciso,klon, klev),     INTENT(OUT)         :: xtu
187  REAL, DIMENSION (ntraciso,klon, klev),     INTENT(OUT)         :: dxtls
188  REAL, DIMENSION (ntraciso,klon, klev),     INTENT(OUT)         :: dxtke
189  REAL, DIMENSION (ntraciso,klon, klev),     INTENT(OUT)         :: d_deltaxtw2
190#endif
191
192  ! Variables internes
193  ! -------------------
194
195  ! Variables à fixer
196  INTEGER, SAVE                                         :: igout
197  !$OMP THREADPRIVATE(igout)
198  LOGICAL, SAVE                                         :: first = .TRUE.
199  !$OMP THREADPRIVATE(first)
200!jyg<
201!!  REAL, SAVE                                            :: stark, wdens_ref, coefgw, alpk
202  REAL, SAVE, DIMENSION(2)                              :: wdens_ref
203  REAL, SAVE                                            :: stark, coefgw, alpk
204!>jyg
205  REAL, SAVE                                            :: crep_upper, crep_sol 
206  !$OMP THREADPRIVATE(stark, wdens_ref, coefgw, alpk, crep_upper, crep_sol)
207
208  REAL, SAVE                                            :: tau_cv
209  !$OMP THREADPRIVATE(tau_cv)
210  REAL, SAVE                                            :: rzero, aa0 ! minimal wake radius and area
211  !$OMP THREADPRIVATE(rzero, aa0)
212
213  LOGICAL, SAVE                                         :: flag_wk_check_trgl
214  !$OMP THREADPRIVATE(flag_wk_check_trgl)
215  INTEGER, SAVE                                         :: iflag_wk_check_trgl
216  !$OMP THREADPRIVATE(iflag_wk_check_trgl)
217  INTEGER, SAVE                                         :: iflag_wk_pop_dyn
218  !$OMP THREADPRIVATE(iflag_wk_pop_dyn)
219
220  REAL                                                  :: delta_t_min
221  INTEGER                                               :: nsub
222  REAL                                                  :: dtimesub
223  REAL, SAVE                                            :: wdensmin
224  !$OMP THREADPRIVATE(wdensmin)
225  REAL, SAVE                                            :: sigmad, hwmin, wapecut, cstart
226  !$OMP THREADPRIVATE(sigmad, hwmin, wapecut, cstart)
227  REAL, SAVE                                            :: sigmaw_max
228  !$OMP THREADPRIVATE(sigmaw_max) 
229  REAL, SAVE                                            :: dens_rate
230  !$OMP THREADPRIVATE(dens_rate)
231  REAL                                                  :: wdens0
232  ! IM 080208
233  LOGICAL, DIMENSION (klon)                             :: gwake
234
235  ! Variables de sauvegarde
236  REAL, DIMENSION (klon, klev)                          :: deltatw0
237  REAL, DIMENSION (klon, klev)                          :: deltaqw0
238  REAL, DIMENSION (klon, klev)                          :: te, qe
239!!  REAL, DIMENSION (klon)                                :: sigmaw1
240#ifdef ISO
241  REAL, DIMENSION (ntraciso,klon, klev)                          :: deltaxtw0
242  REAL, DIMENSION (ntraciso,klon, klev)                          :: xte
243#endif
244
245  ! Variables liees a la dynamique de population
246  REAL, DIMENSION(klon)                                 :: act
247  REAL, DIMENSION(klon)                                 :: rad_wk, tau_wk_inv
248  REAL, DIMENSION(klon)                                 :: f_shear
249  REAL, DIMENSION(klon)                                 :: drdt
250  REAL, DIMENSION(klon)                                 :: d_sig_gen, d_sig_death, d_sig_col
251  REAL, DIMENSION(klon)                                 :: wape1_act, wape2_act
252  LOGICAL, DIMENSION (klon)                             :: kill_wake
253  INTEGER, SAVE                                         :: iflag_wk_act
254  !$OMP THREADPRIVATE(iflag_wk_act)
255  REAL                                                  :: drdt_pos
256  REAL                                                  :: tau_wk_inv_min
257
258  ! Variables pour les GW
259  REAL, DIMENSION (klon)                                :: ll
260  REAL, DIMENSION (klon, klev)                          :: n2
261  REAL, DIMENSION (klon, klev)                          :: cgw
262  REAL, DIMENSION (klon, klev)                          :: tgw
263
264  ! Variables liees au calcul de hw
265  REAL, DIMENSION (klon)                                :: ptop_provis, ptop, ptop_new
266  REAL, DIMENSION (klon)                                :: sum_dth
267  REAL, DIMENSION (klon)                                :: dthmin
268  REAL, DIMENSION (klon)                                :: z, dz, hw0
269  INTEGER, DIMENSION (klon)                             :: ktop, kupper
270
271  ! Variables liees au test de la forme triangulaire du profil de Delta_theta
272  REAL, DIMENSION (klon)                                :: sum_half_dth
273  REAL, DIMENSION (klon)                                :: dz_half
274
275  ! Sub-timestep tendencies and related variables
276  REAL, DIMENSION (klon, klev)                          :: d_deltatw, d_deltaqw
277  REAL, DIMENSION (klon, klev)                          :: d_te, d_qe
278  REAL, DIMENSION (klon)                                :: d_awdens, d_wdens
279  REAL, DIMENSION (klon)                                :: d_sigmaw, alpha
280  REAL, DIMENSION (klon)                                :: q0_min, q1_min
281  LOGICAL, DIMENSION (klon)                             :: wk_adv, ok_qx_qw
282#ifdef ISO
283  REAL, DIMENSION (ntraciso,klon, klev)                          :: d_deltaxtw
284  REAL, DIMENSION (ntraciso,klon, klev)                          :: d_xte
285#endif
286  REAL, SAVE                                            :: epsilon
287  !$OMP THREADPRIVATE(epsilon)
288  DATA epsilon/1.E-15/
289
290  ! Autres variables internes
291  INTEGER                                               ::isubstep, k, i
292
293  REAL                                                  :: wdens_targ
294  REAL                                                  :: sigmaw_targ
295
296  REAL, DIMENSION (klon)                                :: sum_thu, sum_tu, sum_qu, sum_thvu
297  REAL, DIMENSION (klon)                                :: sum_dq, sum_rho
298  REAL, DIMENSION (klon)                                :: sum_dtdwn, sum_dqdwn
299  REAL, DIMENSION (klon)                                :: av_thu, av_tu, av_qu, av_thvu
300  REAL, DIMENSION (klon)                                :: av_dth, av_dq, av_rho
301  REAL, DIMENSION (klon)                                :: av_dtdwn, av_dqdwn
302! pas besoin d'isos ici
303
304  REAL, DIMENSION (klon, klev)                          :: rho, rhow
305  REAL, DIMENSION (klon, klev+1)                        :: rhoh
306  REAL, DIMENSION (klon, klev)                          :: rhow_moyen
307  REAL, DIMENSION (klon, klev)                          :: zh
308  REAL, DIMENSION (klon, klev+1)                        :: zhh
309  REAL, DIMENSION (klon, klev)                          :: epaisseur1, epaisseur2
310
311  REAL, DIMENSION (klon, klev)                          :: the, thu
312
313  REAL, DIMENSION (klon, klev)                          :: omgbw
314  REAL, DIMENSION (klon)                                :: pupper
315  REAL, DIMENSION (klon)                                :: omgtop
316  REAL, DIMENSION (klon, klev)                          :: dp_omgbw
317  REAL, DIMENSION (klon)                                :: ztop, dztop
318  REAL, DIMENSION (klon, klev)                          :: alpha_up
319
320  REAL, DIMENSION (klon)                                :: rre1, rre2
321  REAL                                                  :: rrd1, rrd2
322  REAL, DIMENSION (klon, klev)                          :: th1, th2, q1, q2
323  REAL, DIMENSION (klon, klev)                          :: d_th1, d_th2, d_dth
324  REAL, DIMENSION (klon, klev)                          :: d_q1, d_q2, d_dq
325  REAL, DIMENSION (klon, klev)                          :: omgbdq
326#ifdef ISO
327      REAL, DIMENSION(ntraciso,klon,klev) :: xt1, xt2
328      REAL, DIMENSION(ntraciso,klon,klev) :: D_xt1, D_xt2, D_dxt
329      REAL, DIMENSION(ntraciso,klon,klev) :: omgbdxt
330      integer ixt
331#endif
332
333  REAL, DIMENSION (klon)                                :: ff, gg
334  REAL, DIMENSION (klon)                                :: wape2, cstar2, heff
335                                                       
336  REAL, DIMENSION (klon, klev)                          :: crep
337                                                       
338  REAL, DIMENSION (klon, klev)                          :: ppi
339
340  ! cc nrlmd
341  REAL, DIMENSION (klon)                                :: death_rate
342!!  REAL, DIMENSION (klon)                                :: nat_rate
343  REAL, DIMENSION (klon, klev)                          :: entr
344  REAL, DIMENSION (klon, klev)                          :: detr
345
346  REAL, DIMENSION(klon)                                 :: sigmaw_in             ! pour les prints
347  REAL, DIMENSION(klon)                                 :: awdens_in, wdens_in   ! pour les prints
348
349  ! -------------------------------------------------------------------------
350  ! Initialisations
351  ! -------------------------------------------------------------------------
352
353  ! print*, 'wake initialisations'
354!#ifdef ISOVERIF
355!        write(*,*) 'wake 358: entree'
356!#endif
357
358  ! Essais d'initialisation avec sigmaw = 0.02 et hw = 10.
359  ! -------------------------------------------------------------------------
360
361!!  DATA wapecut, sigmad, hwmin/5., .02, 10./
362!!  DATA wapecut, sigmad, hwmin/1., .02, 10./
363  DATA sigmad, hwmin/.02, 10./
364!!  DATA wdensmin/1.e-12/
365  DATA wdensmin/1.e-14/
366  ! cc nrlmd
367  DATA sigmaw_max/0.4/
368  DATA dens_rate/0.1/
369  ! cc
370  ! Longueur de maille (en m)
371  ! -------------------------------------------------------------------------
372
373  ! ALON = 3.e5
374  ! alon = 1.E6
375
376  !  Provisionnal; to be suppressed when f_shear is parameterized
377  f_shear(:) = 1.       ! 0. for strong shear, 1. for weak shear
378
379
380  ! Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei)
381
382  ! coefgw : Coefficient pour les ondes de gravité
383  ! stark : Coefficient k dans Cstar=k*sqrt(2*WAPE)
384  ! wdens : Densité surfacique de poche froide
385  ! -------------------------------------------------------------------------
386
387  ! cc nrlmd      coefgw=10
388  ! coefgw=1
389  ! wdens0 = 1.0/(alon**2)
390  ! cc nrlmd      wdens = 1.0/(alon**2)
391  ! cc nrlmd      stark = 0.50
392  ! CRtest
393  ! cc nrlmd      alpk=0.1
394  ! alpk = 1.0
395  ! alpk = 0.5
396  ! alpk = 0.05
397
398 if (first) then
399
400  igout = klon/2+1/klon
401
402  crep_upper = 0.9
403  crep_sol = 1.0
404
405  ! Get wapecut from parameter file
406  wapecut = 1.
407  CALL getin_p('wapecut', wapecut)
408
409  ! cc nrlmd Lecture du fichier wake_param.data
410  stark=0.33
411  CALL getin_p('stark',stark)
412  cstart = stark*sqrt(2.*wapecut)
413
414  alpk=0.25
415  CALL getin_p('alpk',alpk)
416!jyg<
417!!  wdens_ref=8.E-12
418!!  CALL getin_p('wdens_ref',wdens_ref)
419  wdens_ref(1)=8.E-12
420  wdens_ref(2)=8.E-12
421  CALL getin_p('wdens_ref_o',wdens_ref(1))    !wake number per unit area ; ocean
422  CALL getin_p('wdens_ref_l',wdens_ref(2))    !wake number per unit area ; land
423!>jyg
424!
425!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
426!!!!!!!!!  Population dynamics parameters    !!!!!!!!!!!!!!!!!!!!!!!!!!!!
427!------------------------------------------------------------------------
428
429  iflag_wk_pop_dyn = 0
430  CALL getin_p('iflag_wk_pop_dyn',iflag_wk_pop_dyn) ! switch between wdens prescribed
431                                                    ! and wdens prognostic
432  iflag_wk_act = 0
433  CALL getin_p('iflag_wk_act',iflag_wk_act) ! 0: act(:)=0.
434                                            ! 1: act(:)=1.
435                                            ! 2: act(:)=f(Wape)
436
437  rzero = 5000.
438  CALL getin_p('rzero_wk', rzero)
439  aa0 = 3.14*rzero*rzero
440!
441  tau_cv = 4000.
442  CALL getin_p('tau_cv', tau_cv)
443
444!------------------------------------------------------------------------
445
446  coefgw=4.
447  CALL getin_p('coefgw',coefgw)
448
449  WRITE(*,*) 'stark=', stark
450  WRITE(*,*) 'alpk=', alpk
451!jyg<
452!!  WRITE(*,*) 'wdens_ref=', wdens_ref
453  WRITE(*,*) 'wdens_ref_o=', wdens_ref(1)
454  WRITE(*,*) 'wdens_ref_l=', wdens_ref(2)
455!>jyg
456  WRITE(*,*) 'iflag_wk_pop_dyn=',iflag_wk_pop_dyn
457  WRITE(*,*) 'iflag_wk_act',iflag_wk_act
458  WRITE(*,*) 'coefgw=', coefgw
459
460  flag_wk_check_trgl=.false.
461  CALL getin_p('flag_wk_check_trgl ', flag_wk_check_trgl)
462  WRITE(*,*) 'flag_wk_check_trgl=', flag_wk_check_trgl
463  WRITE(*,*) 'flag_wk_check_trgl OBSOLETE. Utilisr iflag_wk_check_trgl plutot'
464  iflag_wk_check_trgl=0 ; IF (flag_wk_check_trgl) iflag_wk_check_trgl=1
465  CALL getin_p('iflag_wk_check_trgl ', iflag_wk_check_trgl)
466  WRITE(*,*) 'iflag_wk_check_trgl=', iflag_wk_check_trgl
467
468  first=.false.
469 endif
470
471 IF (iflag_wk_pop_dyn == 0) THEN
472  ! Initialisation de toutes des densites a wdens_ref.
473  ! Les densites peuvent evoluer si les poches debordent
474  ! (voir au tout debut de la boucle sur les substeps)
475  !jyg<
476  !!  wdens(:) = wdens_ref
477   DO i = 1,klon
478     wdens(i) = wdens_ref(znatsurf(i)+1)
479   ENDDO
480  !>jyg
481 ENDIF  ! (iflag_wk_pop_dyn == 0)
482
483  ! print*,'stark',stark
484  ! print*,'alpk',alpk
485  ! print*,'wdens',wdens
486  ! print*,'coefgw',coefgw
487  ! cc
488  ! Minimum value for |T_wake - T_undist|. Used for wake top definition
489  ! -------------------------------------------------------------------------
490
491  delta_t_min = 0.2
492
493  ! 1. - Save initial values, initialize tendencies, initialize output fields
494  ! ------------------------------------------------------------------------
495
496!jyg<
497!!  DO k = 1, klev
498!!    DO i = 1, klon
499!!      ppi(i, k) = pi(i, k)
500!!      deltatw0(i, k) = deltatw(i, k)
501!!      deltaqw0(i, k) = deltaqw(i, k)
502!!      te(i, k) = te0(i, k)
503!!      qe(i, k) = qe0(i, k)
504!!      dtls(i, k) = 0.
505!!      dqls(i, k) = 0.
506!!      d_deltat_gw(i, k) = 0.
507!!      d_te(i, k) = 0.
508!!      d_qe(i, k) = 0.
509!!      d_deltatw(i, k) = 0.
510!!      d_deltaqw(i, k) = 0.
511!!      ! IM 060508 beg
512!!      d_deltatw2(i, k) = 0.
513!!      d_deltaqw2(i, k) = 0.
514!!      ! IM 060508 end
515!!    END DO
516!!  END DO
517      ppi(:,:) = pi(:,:)
518      deltatw0(:,:) = deltatw(:,:)
519      deltaqw0(:,:) = deltaqw(:,:)
520      te(:,:) = te0(:,:)
521      qe(:,:) = qe0(:,:)
522      dtls(:,:) = 0.
523      dqls(:,:) = 0.
524      d_deltat_gw(:,:) = 0.
525      d_te(:,:) = 0.
526      d_qe(:,:) = 0.
527      d_deltatw(:,:) = 0.
528      d_deltaqw(:,:) = 0.
529      d_deltatw2(:,:) = 0.
530      d_deltaqw2(:,:) = 0.
531#ifdef ISO
532      deltaxtw0(:,:,:)= deltaxtw(:,:,:)
533      xte(:,:,:) = xte0(:,:,:)
534      dxtls(:,:,:) = 0.
535      d_xte(:,:,:) = 0.
536      d_deltaxtw(:,:,:) = 0.
537      d_deltaxtw2(:,:,:)=0.
538      xt1(:,:,:) = 0.
539      xt2(:,:,:)=0.   
540      ! init non indispensable mais facilite verif iso
541      q1(:,:)=0.0
542      q2(:,:)=0.0
543#endif
544
545      IF (iflag_wk_act == 0) THEN
546        act(:) = 0.
547      ELSEIF (iflag_wk_act == 1) THEN
548        act(:) = 1.
549      ENDIF
550
551!!  DO i = 1, klon
552!!   sigmaw_in(i) = sigmaw(i)
553!!  END DO
554   sigmaw_in(:) = sigmaw(:)
555!>jyg
556
557  ! sigmaw1=sigmaw
558  ! IF (sigd_con.GT.sigmaw1) THEN
559  ! print*, 'sigmaw,sigd_con', sigmaw, sigd_con
560  ! ENDIF
561  IF (iflag_wk_pop_dyn >=1) THEN
562    DO i = 1, klon
563      wdens_targ = max(wdens(i),wdensmin)
564      d_wdens2(i) = wdens_targ - wdens(i)
565      wdens(i) = wdens_targ
566    END DO
567  ELSE
568    DO i = 1, klon
569      d_awdens2(i) = 0.
570      d_wdens2(i) = 0.
571    END DO
572  ENDIF  ! (iflag_wk_pop_dyn >=1)
573!
574  DO i = 1, klon
575    ! c      sigmaw(i) = amax1(sigmaw(i),sigd_con(i))
576!jyg<
577!!    sigmaw(i) = amax1(sigmaw(i), sigmad)
578!!    sigmaw(i) = amin1(sigmaw(i), 0.99)
579    sigmaw_targ = min(max(sigmaw(i), sigmad),0.99)
580    d_sigmaw2(i) = sigmaw_targ - sigmaw(i)
581    sigmaw(i) = sigmaw_targ
582!>jyg
583  END DO
584
585!
586  IF (iflag_wk_pop_dyn >= 1) THEN
587    awdens_in(:) = awdens(:)
588    wdens_in(:) = wdens(:)
589!!    wdens(:) = wdens(:) + wgen(:)*dtime
590!!    d_wdens2(:) = wgen(:)*dtime
591!!  ELSE
592  ENDIF  ! (iflag_wk_pop_dyn >= 1)
593
594  wape(:) = 0.
595  wape2(:) = 0.
596  d_sigmaw(:) = 0.
597  ktopw(:) = 0
598!
599!<jyg
600dth(:,:) = 0.
601tu(:,:) = 0.
602qu(:,:) = 0.
603dtke(:,:) = 0.
604dqke(:,:) = 0.
605spread(:,:) = 0.
606omgbdth(:,:) = 0.
607omg(:,:) = 0.
608dp_omgb(:,:) = 0.
609dp_deltomg(:,:) = 0.
610hw(:) = 0.
611wape(:) = 0.
612fip(:) = 0.
613gfl(:) = 0.
614cstar(:) = 0.
615ktopw(:) = 0
616!
617!  Vertical advection local variables
618omgbw(:,:) = 0.
619omgtop(:) = 0
620dp_omgbw(:,:) = 0.
621omgbdq(:,:) = 0.
622#ifdef ISO
623xtu(:,:,:) = 0.
624dxtke(:,:,:) = 0.
625omgbdxt(:,:,:) = 0.
626#endif
627!>jyg
628!
629  IF (prt_level>=10) THEN
630    PRINT *, 'wake-1, sigmaw(igout) ', sigmaw(igout)
631    PRINT *, 'wake-1, deltatw(igout,k) ', (k,deltatw(igout,k), k=1,klev)
632    PRINT *, 'wake-1, deltaqw(igout,k) ', (k,deltaqw(igout,k), k=1,klev)
633    PRINT *, 'wake-1, dowwdraughts, amdwn(igout,k) ', (k,amdwn(igout,k), k=1,klev)
634    PRINT *, 'wake-1, dowwdraughts, dtdwn(igout,k) ', (k,dtdwn(igout,k), k=1,klev)
635    PRINT *, 'wake-1, dowwdraughts, dqdwn(igout,k) ', (k,dqdwn(igout,k), k=1,klev)
636    PRINT *, 'wake-1, updraughts, amup(igout,k) ', (k,amup(igout,k), k=1,klev)
637    PRINT *, 'wake-1, updraughts, dta(igout,k) ', (k,dta(igout,k), k=1,klev)
638    PRINT *, 'wake-1, updraughts, dqa(igout,k) ', (k,dqa(igout,k), k=1,klev)
639  ENDIF
640
641  ! 2. - Prognostic part
642  ! --------------------
643
644
645  ! 2.1 - Undisturbed area and Wake integrals
646  ! ---------------------------------------------------------
647
648  DO i = 1, klon
649    z(i) = 0.
650    ktop(i) = 0
651    kupper(i) = 0
652    sum_thu(i) = 0.
653    sum_tu(i) = 0.
654    sum_qu(i) = 0.
655    sum_thvu(i) = 0.
656    sum_dth(i) = 0.
657    sum_dq(i) = 0.
658    sum_rho(i) = 0.
659    sum_dtdwn(i) = 0.
660    sum_dqdwn(i) = 0.
661
662    av_thu(i) = 0.
663    av_tu(i) = 0.
664    av_qu(i) = 0.
665    av_thvu(i) = 0.
666    av_dth(i) = 0.
667    av_dq(i) = 0.
668    av_rho(i) = 0.
669    av_dtdwn(i) = 0.
670    av_dqdwn(i) = 0.
671    ! pas besoin d'isos ici
672  END DO
673
674
675#ifdef ISOVERIF
676      ! TODO     
677#endif
678
679  ! Distance between wakes
680  DO i = 1, klon
681    ll(i) = (1-sqrt(sigmaw(i)))/sqrt(wdens(i))
682  END DO
683  ! Potential temperatures and humidity
684  ! ----------------------------------------------------------
685  DO k = 1, klev
686    DO i = 1, klon
687      ! write(*,*)'wake 1',i,k,rd,te(i,k)
688      rho(i, k) = p(i, k)/(rd*te(i,k))
689      ! write(*,*)'wake 2',rho(i,k)
690      IF (k==1) THEN
691        ! write(*,*)'wake 3',i,k,rd,te(i,k)
692        rhoh(i, k) = ph(i, k)/(rd*te(i,k))
693        ! write(*,*)'wake 4',i,k,rd,te(i,k)
694        zhh(i, k) = 0
695      ELSE
696        ! write(*,*)'wake 5',rd,(te(i,k)+te(i,k-1))
697        rhoh(i, k) = ph(i, k)*2./(rd*(te(i,k)+te(i,k-1)))
698        ! write(*,*)'wake 6',(-rhoh(i,k)*RG)+zhh(i,k-1)
699        zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*rg) + zhh(i, k-1)
700      END IF
701      ! write(*,*)'wake 7',ppi(i,k)
702      the(i, k) = te(i, k)/ppi(i, k)
703      thu(i, k) = (te(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k)
704      tu(i, k) = te(i, k) - deltatw(i, k)*sigmaw(i)
705      qu(i, k) = qe(i, k) - deltaqw(i, k)*sigmaw(i)
706      ! write(*,*)'wake 8',(rd*(te(i,k)+deltatw(i,k)))
707      rhow(i, k) = p(i, k)/(rd*(te(i,k)+deltatw(i,k)))
708      dth(i, k) = deltatw(i, k)/ppi(i, k)
709#ifdef ISO
710        do ixt=1,ntraciso
711          xtu(ixt,i,k)  =  xte(ixt,i,k) - deltaxtw(ixt,i,k)*sigmaw(i)
712        enddo !do ixt=1,ntraciso
713#ifdef ISOVERIF
714        if (iso_eau.gt.0) then
715              call iso_verif_egalite(deltaqw(i,k),deltaxtw(iso_eau,i,k),'wake 723a')
716              call iso_verif_egalite(qu(i,k),xtu(iso_eau,i,k),'wake 723b')
717        endif
718        if (iso_HDO.gt.0) then
719               if (iso_verif_aberrant_enc_choix_nostop(xtu(iso_hdo,i,k),qu(i,k),ridicule,deltalim, &
720                        'wake 723c xtu').eq.1) then
721                     stop
722               endif
723         endif 
724#endif
725#endif
726    END DO
727  END DO
728
729
730
731  DO k = 1, klev - 1
732    DO i = 1, klon
733      IF (k==1) THEN
734        n2(i, k) = 0
735      ELSE
736        n2(i, k) = amax1(0., -rg**2/the(i,k)*rho(i,k)*(the(i,k+1)-the(i,k-1))/ &
737                             (p(i,k+1)-p(i,k-1)))
738      END IF
739      zh(i, k) = (zhh(i,k)+zhh(i,k+1))/2
740
741      cgw(i, k) = sqrt(n2(i,k))*zh(i, k)
742      tgw(i, k) = coefgw*cgw(i, k)/ll(i)
743    END DO
744  END DO
745
746  DO i = 1, klon
747    n2(i, klev) = 0
748    zh(i, klev) = 0
749    cgw(i, klev) = 0
750    tgw(i, klev) = 0
751  END DO
752
753  ! Calcul de la masse volumique moyenne de la colonne   (bdlmd)
754  ! -----------------------------------------------------------------
755
756  DO k = 1, klev
757    DO i = 1, klon
758      epaisseur1(i, k) = 0.
759      epaisseur2(i, k) = 0.
760    END DO
761  END DO
762
763  DO i = 1, klon
764    epaisseur1(i, 1) = -(ph(i,2)-ph(i,1))/(rho(i,1)*rg) + 1.
765    epaisseur2(i, 1) = -(ph(i,2)-ph(i,1))/(rho(i,1)*rg) + 1.
766    rhow_moyen(i, 1) = rhow(i, 1)
767  END DO
768
769  DO k = 2, klev
770    DO i = 1, klon
771      epaisseur1(i, k) = -(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg) + 1.
772      epaisseur2(i, k) = epaisseur2(i, k-1) + epaisseur1(i, k)
773      rhow_moyen(i, k) = (rhow_moyen(i,k-1)*epaisseur2(i,k-1)+rhow(i,k)* &
774        epaisseur1(i,k))/epaisseur2(i, k)
775    END DO
776  END DO
777
778
779  ! Choose an integration bound well above wake top
780  ! -----------------------------------------------------------------
781
782  ! Pupper = 50000.  ! melting level
783  ! Pupper = 60000.
784  ! Pupper = 80000.  ! essais pour case_e
785  DO i = 1, klon
786    pupper(i) = 0.6*ph(i, 1)
787    pupper(i) = max(pupper(i), 45000.)
788    ! cc        Pupper(i) = 60000.
789  END DO
790
791
792  ! Determine Wake top pressure (Ptop) from buoyancy integral
793  ! --------------------------------------------------------
794
795  ! -1/ Pressure of the level where dth becomes less than delta_t_min.
796
797  DO i = 1, klon
798    ptop_provis(i) = ph(i, 1)
799  END DO
800  DO k = 2, klev
801    DO i = 1, klon
802
803      ! IM v3JYG; ptop_provis(i).LT. ph(i,1)
804
805      IF (dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min .AND. &
806          ptop_provis(i)==ph(i,1)) THEN
807        ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)- &
808                          (dth(i,k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
809      END IF
810    END DO
811  END DO
812
813  ! -2/ dth integral
814
815  DO i = 1, klon
816    sum_dth(i) = 0.
817    dthmin(i) = -delta_t_min
818    z(i) = 0.
819  END DO
820
821  DO k = 1, klev
822    DO i = 1, klon
823      dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*rg)
824      IF (dz(i)>0) THEN
825        z(i) = z(i) + dz(i)
826        sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
827        dthmin(i) = amin1(dthmin(i), dth(i,k))
828      END IF
829    END DO
830  END DO
831
832  ! -3/ height of triangle with area= sum_dth and base = dthmin
833
834  DO i = 1, klon
835    hw0(i) = 2.*sum_dth(i)/amin1(dthmin(i), -0.5)
836    hw0(i) = amax1(hwmin, hw0(i))
837  END DO
838
839  ! -4/ now, get Ptop
840
841  DO i = 1, klon
842    z(i) = 0.
843    ptop(i) = ph(i, 1)
844  END DO
845
846  DO k = 1, klev
847    DO i = 1, klon
848      dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg), hw0(i)-z(i))
849      IF (dz(i)>0) THEN
850        z(i) = z(i) + dz(i)
851        ptop(i) = ph(i, k) - rho(i, k)*rg*dz(i)
852      END IF
853    END DO
854  END DO
855
856  IF (prt_level>=10) THEN
857    PRINT *, 'wake-2, ptop_provis(igout), ptop(igout) ', ptop_provis(igout), ptop(igout)
858  ENDIF
859
860
861  ! -5/ Determination de ktop et kupper
862
863  DO k = klev, 1, -1
864    DO i = 1, klon
865      IF (ph(i,k+1)<ptop(i)) ktop(i) = k
866      IF (ph(i,k+1)<pupper(i)) kupper(i) = k
867    END DO
868  END DO
869
870  ! On evite kupper = 1 et kupper = klev
871  DO i = 1, klon
872    kupper(i) = max(kupper(i), 2)
873    kupper(i) = min(kupper(i), klev-1)
874  END DO
875
876
877  ! -6/ Correct ktop and ptop
878
879  DO i = 1, klon
880    ptop_new(i) = ptop(i)
881  END DO
882  DO k = klev, 2, -1
883    DO i = 1, klon
884      IF (k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. &
885          dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
886        ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, &
887          k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
888      END IF
889    END DO
890  END DO
891
892  DO i = 1, klon
893    ptop(i) = ptop_new(i)
894  END DO
895
896  DO k = klev, 1, -1
897    DO i = 1, klon
898      IF (ph(i,k+1)<ptop(i)) ktop(i) = k
899    END DO
900  END DO
901
902  IF (prt_level>=10) THEN
903    PRINT *, 'wake-3, ktop(igout), kupper(igout) ', ktop(igout), kupper(igout)
904  ENDIF
905
906  ! -5/ Set deltatw & deltaqw to 0 above kupper
907
908  DO k = 1, klev
909    DO i = 1, klon
910      IF (k>=kupper(i)) THEN
911        deltatw(i, k) = 0.
912        deltaqw(i, k) = 0.
913        d_deltatw2(i,k) = -deltatw0(i,k)
914#ifdef ISO
915        do ixt=1,ntraciso
916        deltaxtw(ixt,i, k) = 0.
917        d_deltaxtw2(ixt,i,k) = -deltaxtw0(ixt,i,k)
918        enddo
919#endif
920      END IF
921    END DO
922  END DO
923
924
925  ! Vertical gradient of LS omega
926
927  DO k = 1, klev
928    DO i = 1, klon
929      IF (k<=kupper(i)) THEN
930        dp_omgb(i, k) = (omgb(i,k+1)-omgb(i,k))/(ph(i,k+1)-ph(i,k))
931      END IF
932    END DO
933  END DO
934
935  ! Integrals (and wake top level number)
936  ! --------------------------------------
937
938  ! Initialize sum_thvu to 1st level virt. pot. temp.
939
940  DO i = 1, klon
941    z(i) = 1.
942    dz(i) = 1.
943    sum_thvu(i) = thu(i, 1)*(1.+epsim1*qu(i,1))*dz(i)
944    sum_dth(i) = 0.
945  END DO
946
947  DO k = 1, klev
948    DO i = 1, klon
949      dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
950      IF (dz(i)>0) THEN
951        z(i) = z(i) + dz(i)
952        sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i)
953        sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i)
954        sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i)
955        sum_thvu(i) = sum_thvu(i) + thu(i, k)*(1.+epsim1*qu(i,k))*dz(i)
956        sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
957        sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
958        sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i)
959        sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
960        sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
961      END IF
962    END DO
963  END DO
964
965  DO i = 1, klon
966    hw0(i) = z(i)
967  END DO
968
969
970  ! 2.1 - WAPE and mean forcing computation
971  ! ---------------------------------------
972
973  ! ---------------------------------------
974
975  ! Means
976
977  DO i = 1, klon
978    av_thu(i) = sum_thu(i)/hw0(i)
979    av_tu(i) = sum_tu(i)/hw0(i)
980    av_qu(i) = sum_qu(i)/hw0(i)
981    av_thvu(i) = sum_thvu(i)/hw0(i)
982    ! av_thve = sum_thve/hw0
983    av_dth(i) = sum_dth(i)/hw0(i)
984    av_dq(i) = sum_dq(i)/hw0(i)
985    av_rho(i) = sum_rho(i)/hw0(i)
986    av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
987    av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
988
989    wape(i) = -rg*hw0(i)*(av_dth(i)+ &
990        epsim1*(av_thu(i)*av_dq(i)+av_dth(i)*av_qu(i)+av_dth(i)*av_dq(i)))/av_thvu(i)
991
992  END DO
993
994  ! 2.2 Prognostic variable update
995  ! ------------------------------
996
997  ! Filter out bad wakes
998
999  DO k = 1, klev
1000    DO i = 1, klon
1001      IF (wape(i)<0.) THEN
1002        deltatw(i, k) = 0.
1003        deltaqw(i, k) = 0.
1004        dth(i, k) = 0.
1005        d_deltatw2(i,k) = -deltatw0(i,k)
1006        d_deltaqw2(i,k) = -deltaqw0(i,k)
1007#ifdef ISO
1008        do ixt=1,ntraciso
1009        deltaxtw(ixt,i, k) = 0.
1010        d_deltaxtw2(ixt,i,k) = -deltaxtw0(ixt,i,k)
1011        enddo
1012#endif
1013      END IF
1014    END DO
1015  END DO
1016
1017  DO i = 1, klon
1018    IF (wape(i)<0.) THEN
1019      wape(i) = 0.
1020      cstar(i) = 0.
1021      hw(i) = hwmin
1022!jyg<
1023!!      sigmaw(i) = amax1(sigmad, sigd_con(i))
1024      sigmaw_targ = max(sigmad, sigd_con(i))
1025      d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
1026      sigmaw(i) = sigmaw_targ
1027!>jyg
1028      fip(i) = 0.
1029      gwake(i) = .FALSE.
1030    ELSE
1031      hw(i) = hw0(i)
1032      cstar(i) = stark*sqrt(2.*wape(i))
1033      gwake(i) = .TRUE.
1034    END IF
1035  END DO
1036
1037
1038  ! Check qx and qw positivity
1039  ! --------------------------
1040  DO i = 1, klon
1041    q0_min(i) = min((qe(i,1)-sigmaw(i)*deltaqw(i,1)),  &
1042                    (qe(i,1)+(1.-sigmaw(i))*deltaqw(i,1)))
1043  END DO
1044  DO k = 2, klev
1045    DO i = 1, klon
1046      q1_min(i) = min((qe(i,k)-sigmaw(i)*deltaqw(i,k)), &
1047                      (qe(i,k)+(1.-sigmaw(i))*deltaqw(i,k)))
1048      IF (q1_min(i)<=q0_min(i)) THEN
1049        q0_min(i) = q1_min(i)
1050      END IF
1051    END DO
1052  END DO
1053
1054  DO i = 1, klon
1055    ok_qx_qw(i) = q0_min(i) >= 0.
1056    alpha(i) = 1.
1057  END DO
1058
1059  IF (prt_level>=10) THEN
1060    PRINT *, 'wake-4, sigmaw(igout), cstar(igout), wape(igout), ktop(igout) ', &
1061                      sigmaw(igout), cstar(igout), wape(igout), ktop(igout)
1062  ENDIF
1063
1064
1065  ! C -----------------------------------------------------------------
1066  ! Sub-time-stepping
1067  ! -----------------
1068
1069  nsub = 10
1070  dtimesub = dtime/nsub
1071
1072  ! ------------------------------------------------------------
1073  DO isubstep = 1, nsub
1074    ! ------------------------------------------------------------
1075
1076    ! wk_adv is the logical flag enabling wake evolution in the time advance
1077    ! loop
1078    DO i = 1, klon
1079      wk_adv(i) = ok_qx_qw(i) .AND. alpha(i) >= 1.
1080    END DO
1081    IF (prt_level>=10) THEN
1082      PRINT *, 'wake-4.1, isubstep,wk_adv(igout),cstar(igout),wape(igout), ptop(igout) ', &
1083                          isubstep,wk_adv(igout),cstar(igout),wape(igout), ptop(igout)
1084    ENDIF
1085
1086    ! cc nrlmd   Ajout d'un recalcul de wdens dans le cas d'un entrainement
1087    ! négatif de ktop à kupper --------
1088    ! cc           On calcule pour cela une densité wdens0 pour laquelle on
1089    ! aurait un entrainement nul ---
1090    !jyg<
1091    ! Dans la configuration avec wdens prognostique, il s'agit d'un cas ou
1092    ! les poches sont insuffisantes pour accueillir tout le flux de masse
1093    ! des descentes unsaturees. Nous faisons alors l'hypothese que la
1094    ! convection profonde cree directement de nouvelles poches, sans passer
1095    ! par les thermiques. La nouvelle valeur de wdens est alors imposée.
1096
1097    DO i = 1, klon
1098      ! c       print *,' isubstep,wk_adv(i),cstar(i),wape(i) ',
1099      ! c     $           isubstep,wk_adv(i),cstar(i),wape(i)
1100      IF (wk_adv(i) .AND. cstar(i)>0.01) THEN
1101        omg(i, kupper(i)+1) = -rg*amdwn(i, kupper(i)+1)/sigmaw(i) + &
1102                               rg*amup(i, kupper(i)+1)/(1.-sigmaw(i))
1103        wdens0 = (sigmaw(i)/(4.*3.14))* &
1104          ((1.-sigmaw(i))*omg(i,kupper(i)+1)/((ph(i,1)-pupper(i))*cstar(i)))**(2)
1105        IF (prt_level >= 10) THEN
1106             print*,'omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i), ph(i,1)-pupper(i)', &
1107                     omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i), ph(i,1)-pupper(i)
1108        ENDIF
1109        IF (wdens(i)<=wdens0*1.1) THEN
1110          IF (iflag_wk_pop_dyn >= 1) THEN
1111             d_wdens2(i) = d_wdens2(i) + wdens0 - wdens(i)
1112          ENDIF
1113          wdens(i) = wdens0
1114        END IF
1115      END IF
1116    END DO
1117
1118    DO i = 1, klon
1119      IF (wk_adv(i)) THEN
1120        gfl(i) = 2.*sqrt(3.14*wdens(i)*sigmaw(i))
1121        rad_wk(i) = sqrt(sigmaw(i)/(3.14*wdens(i)))
1122!jyg<
1123!!        sigmaw(i) = amin1(sigmaw(i), sigmaw_max)
1124        sigmaw_targ = min(sigmaw(i), sigmaw_max)
1125        d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
1126        sigmaw(i) = sigmaw_targ
1127!>jyg
1128      END IF
1129    END DO
1130
1131    IF (iflag_wk_pop_dyn >= 1) THEN
1132!    The variable "death_rate" is significant only when iflag_wk_pop_dyn = 0.
1133!    Here, it has to be set to zero.
1134      death_rate(:) = 0.
1135
1136      IF (iflag_wk_act ==2) THEN
1137      DO i = 1, klon
1138        IF (wk_adv(i)) THEN
1139          wape1_act(i) = abs(cin(i))
1140          wape2_act(i) = 2.*wape1_act(i) + 1.
1141          act(i) = min(1., max(0., (wape(i)-wape1_act(i)) / (wape2_act(i)-wape1_act(i)) ))
1142        ENDIF  ! (wk_adv(i))
1143      ENDDO
1144      ENDIF  ! (iflag_wk_act ==2)
1145
1146
1147      DO i = 1, klon
1148        IF (wk_adv(i)) THEN
1149!!          tau_wk(i) = max(rad_wk(i)/(3.*cstar(i))*((cstar(i)/cstart)**1.5 - 1), 100.)
1150          tau_wk_inv(i) = max( (3.*cstar(i))/(rad_wk(i)*((cstar(i)/cstart)**1.5 - 1)), 0.)
1151          tau_wk_inv_min = min(tau_wk_inv(i), 1./dtimesub)
1152          drdt(i) = (cstar(i) - wgen(i)*(sigmaw(i)/wdens(i)-aa0)/gfl(i)) / &
1153                    (1 + 2*f_shear(i)*(2.*sigmaw(i)-aa0*wdens(i)) - 2.*sigmaw(i))
1154!!                    (1 - 2*sigmaw(i)*(1.-f_shear(i)))
1155          drdt_pos=max(drdt(i),0.)
1156
1157!!          d_wdens(i) = ( wgen(i)*(1.+2.*(sigmaw(i)-sigmad)) &
1158!!                     - wdens(i)*tau_wk_inv_min &
1159!!                     - 2.*gfl(i)*wdens(i)*Cstar(i) )*dtimesub
1160          d_awdens(i) = ( wgen(i) - (1./tau_cv)*(awdens(i) - act(i)*wdens(i)) )*dtimesub
1161          d_wdens(i) = ( wgen(i) - (wdens(i)-awdens(i))*tau_wk_inv_min -  &
1162                         2.*wdens(i)*gfl(i)*drdt_pos )*dtimesub
1163          d_wdens(i) = max(d_wdens(i), wdensmin-wdens(i))
1164
1165!!          d_sigmaw(i) = ( (1.-2*f_shear(i)*sigmaw(i))*(gfl(i)*Cstar(i)+wgen(i)*sigmad/wdens(i)) &
1166!!                      + 2.*f_shear(i)*wgen(i)*sigmaw(i)**2/wdens(i) &
1167!!                      - sigmaw(i)*tau_wk_inv_min )*dtimesub
1168          d_sig_gen(i) = wgen(i)*aa0
1169          d_sig_death(i) = - sigmaw(i)*(1.-awdens(i)/wdens(i))*tau_wk_inv_min
1170!!          d_sig_col(i) = - 2*f_shear(i)*sigmaw(i)*gfl(i)*drdt_pos
1171          d_sig_col(i) = - 2*f_shear(i)*(2.*sigmaw(i)-wdens(i)*aa0)*gfl(i)*drdt_pos
1172          d_sigmaw(i) = ( d_sig_gen(i) + d_sig_death(i) + d_sig_col(i) + gfl(i)*cstar(i) )*dtimesub
1173          d_sigmaw(i) = max(d_sigmaw(i), sigmad-sigmaw(i))
1174        ENDIF
1175      ENDDO
1176
1177      IF (prt_level >= 10) THEN
1178        print *,'wake, cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), drdt(1) ', &
1179                       cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), drdt(1)
1180        print *,'wake, wdens(1), awdens(1), act(1), d_awdens(1) ', &
1181                       wdens(1), awdens(1), act(1), d_awdens(1)
1182        print *,'wake, wgen, -(wdens-awdens)*tau_wk_inv, -2.*wdens*gfl*drdt_pos, d_wdens ', &
1183                       wgen(1), -(wdens(1)-awdens(1))*tau_wk_inv(1), -2.*wdens(1)*gfl(1)*drdt_pos, d_wdens(1)
1184        print *,'wake, d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1) ', &
1185                       d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1)
1186      ENDIF
1187   
1188    ELSE  !  (iflag_wk_pop_dyn >= 1)
1189
1190    ! cc nrlmd
1191
1192      DO i = 1, klon
1193        IF (wk_adv(i)) THEN
1194          ! cc nrlmd          Introduction du taux de mortalité des poches et
1195          ! test sur sigmaw_max=0.4
1196          ! cc         d_sigmaw(i) = gfl(i)*Cstar(i)*dtimesub
1197          IF (sigmaw(i)>=sigmaw_max) THEN
1198            death_rate(i) = gfl(i)*cstar(i)/sigmaw(i)
1199          ELSE
1200            death_rate(i) = 0.
1201          END IF
1202   
1203          d_sigmaw(i) = gfl(i)*cstar(i)*dtimesub - death_rate(i)*sigmaw(i)* &
1204            dtimesub
1205          ! $              - nat_rate(i)*sigmaw(i)*dtimesub
1206          ! c        print*, 'd_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
1207          ! c     $  death_rate(i),ktop(i),kupper(i)',
1208          ! c     $              d_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
1209          ! c     $  death_rate(i),ktop(i),kupper(i)
1210   
1211          ! sigmaw(i) =sigmaw(i) + gfl(i)*Cstar(i)*dtimesub
1212          ! sigmaw(i) =min(sigmaw(i),0.99)     !!!!!!!!
1213          ! wdens = wdens0/(10.*sigmaw)
1214          ! sigmaw =max(sigmaw,sigd_con)
1215          ! sigmaw =max(sigmaw,sigmad)
1216        END IF
1217      END DO
1218
1219    ENDIF   !  (iflag_wk_pop_dyn >= 1)
1220
1221
1222    ! calcul de la difference de vitesse verticale poche - zone non perturbee
1223    ! IM 060208 differences par rapport au code initial; init. a 0 dp_deltomg
1224    ! IM 060208 et omg sur les niveaux de 1 a klev+1, alors que avant l'on definit
1225    ! IM 060208 au niveau k=1..?
1226    !JYG 161013 Correction : maintenant omg est dimensionne a klev.
1227    DO k = 1, klev
1228      DO i = 1, klon
1229        IF (wk_adv(i)) THEN !!! nrlmd
1230          dp_deltomg(i, k) = 0.
1231        END IF
1232      END DO
1233    END DO
1234    DO k = 1, klev
1235      DO i = 1, klon
1236        IF (wk_adv(i)) THEN !!! nrlmd
1237          omg(i, k) = 0.
1238        END IF
1239      END DO
1240    END DO
1241
1242    DO i = 1, klon
1243      IF (wk_adv(i)) THEN
1244        z(i) = 0.
1245        omg(i, 1) = 0.
1246        dp_deltomg(i, 1) = -(gfl(i)*cstar(i))/(sigmaw(i)*(1-sigmaw(i)))
1247      END IF
1248    END DO
1249
1250    DO k = 2, klev
1251      DO i = 1, klon
1252        IF (wk_adv(i) .AND. k<=ktop(i)) THEN
1253          dz(i) = -(ph(i,k)-ph(i,k-1))/(rho(i,k-1)*rg)
1254          z(i) = z(i) + dz(i)
1255          dp_deltomg(i, k) = dp_deltomg(i, 1)
1256          omg(i, k) = dp_deltomg(i, 1)*z(i)
1257        END IF
1258      END DO
1259    END DO
1260
1261    DO i = 1, klon
1262      IF (wk_adv(i)) THEN
1263        dztop(i) = -(ptop(i)-ph(i,ktop(i)))/(rho(i,ktop(i))*rg)
1264        ztop(i) = z(i) + dztop(i)
1265        omgtop(i) = dp_deltomg(i, 1)*ztop(i)
1266      END IF
1267    END DO
1268
1269    IF (prt_level>=10) THEN
1270      PRINT *, 'wake-4.2, omg(igout,k) ', (k,omg(igout,k), k=1,klev)
1271      PRINT *, 'wake-4.2, omgtop(igout), ptop(igout), ktop(igout) ', &
1272                          omgtop(igout), ptop(igout), ktop(igout)
1273    ENDIF
1274
1275    ! -----------------
1276    ! From m/s to Pa/s
1277    ! -----------------
1278
1279    DO i = 1, klon
1280      IF (wk_adv(i)) THEN
1281        omgtop(i) = -rho(i, ktop(i))*rg*omgtop(i)
1282        dp_deltomg(i, 1) = omgtop(i)/(ptop(i)-ph(i,1))
1283      END IF
1284    END DO
1285
1286    DO k = 1, klev
1287      DO i = 1, klon
1288        IF (wk_adv(i) .AND. k<=ktop(i)) THEN
1289          omg(i, k) = -rho(i, k)*rg*omg(i, k)
1290          dp_deltomg(i, k) = dp_deltomg(i, 1)
1291        END IF
1292      END DO
1293    END DO
1294
1295    ! raccordement lineaire de omg de ptop a pupper
1296
1297    DO i = 1, klon
1298      IF (wk_adv(i) .AND. kupper(i)>ktop(i)) THEN
1299        omg(i, kupper(i)+1) = -rg*amdwn(i, kupper(i)+1)/sigmaw(i) + &
1300          rg*amup(i, kupper(i)+1)/(1.-sigmaw(i))
1301        dp_deltomg(i, kupper(i)) = (omgtop(i)-omg(i,kupper(i)+1))/ &
1302          (ptop(i)-pupper(i))
1303      END IF
1304    END DO
1305
1306    ! c      DO i=1,klon
1307    ! c        print*,'Pente entre 0 et kupper (référence)'
1308    ! c     $           ,omg(i,kupper(i)+1)/(pupper(i)-ph(i,1))
1309    ! c        print*,'Pente entre ktop et kupper'
1310    ! c     $   ,(omg(i,kupper(i)+1)-omgtop(i))/(pupper(i)-ptop(i))
1311    ! c      ENDDO
1312    ! c
1313    DO k = 1, klev
1314      DO i = 1, klon
1315        IF (wk_adv(i) .AND. k>ktop(i) .AND. k<=kupper(i)) THEN
1316          dp_deltomg(i, k) = dp_deltomg(i, kupper(i))
1317          omg(i, k) = omgtop(i) + (ph(i,k)-ptop(i))*dp_deltomg(i, kupper(i))
1318        END IF
1319      END DO
1320    END DO
1321!!    print *,'omg(igout,k) ', (k,omg(igout,k),k=1,klev)
1322    ! cc nrlmd
1323    ! c      DO i=1,klon
1324    ! c      print*,'deltaw_ktop,deltaw_conv',omgtop(i),omg(i,kupper(i)+1)
1325    ! c      END DO
1326    ! cc
1327
1328
1329    ! --    Compute wake average vertical velocity omgbw
1330
1331
1332    DO k = 1, klev
1333      DO i = 1, klon
1334        IF (wk_adv(i)) THEN
1335          omgbw(i, k) = omgb(i, k) + (1.-sigmaw(i))*omg(i, k)
1336        END IF
1337      END DO
1338    END DO
1339    ! --    and its vertical gradient dp_omgbw
1340
1341    DO k = 1, klev-1
1342      DO i = 1, klon
1343        IF (wk_adv(i)) THEN
1344          dp_omgbw(i, k) = (omgbw(i,k+1)-omgbw(i,k))/(ph(i,k+1)-ph(i,k))
1345        END IF
1346      END DO
1347    END DO
1348    DO i = 1, klon
1349      IF (wk_adv(i)) THEN
1350          dp_omgbw(i, klev) = 0.
1351      END IF
1352    END DO
1353
1354    ! --    Upstream coefficients for omgb velocity
1355    ! --    (alpha_up(k) is the coefficient of the value at level k)
1356    ! --    (1-alpha_up(k) is the coefficient of the value at level k-1)
1357    DO k = 1, klev
1358      DO i = 1, klon
1359        IF (wk_adv(i)) THEN
1360          alpha_up(i, k) = 0.
1361          IF (omgb(i,k)>0.) alpha_up(i, k) = 1.
1362        END IF
1363      END DO
1364    END DO
1365
1366    ! Matrix expressing [The,deltatw] from  [Th1,Th2]
1367
1368    DO i = 1, klon
1369      IF (wk_adv(i)) THEN
1370        rre1(i) = 1. - sigmaw(i)
1371        rre2(i) = sigmaw(i)
1372      END IF
1373    END DO
1374    rrd1 = -1.
1375    rrd2 = 1.
1376
1377    ! --    Get [Th1,Th2], dth and [q1,q2]
1378
1379    DO k = 1, klev
1380      DO i = 1, klon
1381        IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN
1382          dth(i, k) = deltatw(i, k)/ppi(i, k)
1383          th1(i, k) = the(i, k) - sigmaw(i)*dth(i, k) ! undisturbed area
1384          th2(i, k) = the(i, k) + (1.-sigmaw(i))*dth(i, k) ! wake
1385          q1(i, k) = qe(i, k) - sigmaw(i)*deltaqw(i, k) ! undisturbed area
1386          q2(i, k) = qe(i, k) + (1.-sigmaw(i))*deltaqw(i, k) ! wake
1387#ifdef ISO
1388          do ixt=1,ntraciso
1389            xt1(ixt,i,k) = xte(ixt,i,k) -sigmaw(i)     *deltaxtw(ixt,i,k) ! undisturbed area
1390            xt2(ixt,i,k) = xte(ixt,i,k) +(1.-sigmaw(i))*deltaxtw(ixt,i,k) ! wake
1391          enddo
1392#endif
1393        END IF
1394      END DO
1395    END DO
1396
1397    DO i = 1, klon
1398      IF (wk_adv(i)) THEN !!! nrlmd
1399        d_th1(i, 1) = 0.
1400        d_th2(i, 1) = 0.
1401        d_dth(i, 1) = 0.
1402        d_q1(i, 1) = 0.
1403        d_q2(i, 1) = 0.
1404        d_dq(i, 1) = 0.
1405#ifdef ISO
1406        do ixt=1,ntraciso
1407          d_xt1(ixt,i,1) = 0.
1408          d_xt2(ixt,i,1) = 0.
1409          d_dxt(ixt,i,1) = 0.
1410        enddo !do ixt=1,ntraciso
1411#endif
1412      END IF
1413    END DO
1414
1415    DO k = 2, klev
1416      DO i = 1, klon
1417        IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN
1418          d_th1(i, k) = th1(i, k-1) - th1(i, k)
1419          d_th2(i, k) = th2(i, k-1) - th2(i, k)
1420          d_dth(i, k) = dth(i, k-1) - dth(i, k)
1421          d_q1(i, k) = q1(i, k-1) - q1(i, k)
1422          d_q2(i, k) = q2(i, k-1) - q2(i, k)
1423          d_dq(i, k) = deltaqw(i, k-1) - deltaqw(i, k)
1424#ifdef ISO
1425        do ixt=1,ntraciso
1426          d_xt1(ixt,i,k) = xt1(ixt,i,k-1)-xt1(ixt,i,k)
1427          d_xt2(ixt,i,k) = xt2(ixt,i,k-1)-xt2(ixt,i,k)
1428          d_dxt(ixt,i,k) = deltaxtw(ixt,i,k-1)-deltaxtw(ixt,i,k)
1429        enddo !do ixt=1,ntraciso
1430#endif
1431        END IF
1432      END DO
1433    END DO
1434
1435    DO i = 1, klon
1436      IF (wk_adv(i)) THEN
1437        omgbdth(i, 1) = 0.
1438        omgbdq(i, 1) = 0.
1439#ifdef ISO
1440        do ixt=1,ntraciso
1441         omgbdxt(ixt,i,1) = 0.
1442        enddo !do ixt=1,ntraciso
1443#endif
1444      END IF
1445    END DO
1446
1447    DO k = 2, klev
1448      DO i = 1, klon
1449        IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN !   loop on interfaces
1450          omgbdth(i, k) = omgb(i, k)*(dth(i,k-1)-dth(i,k))
1451          omgbdq(i, k) = omgb(i, k)*(deltaqw(i,k-1)-deltaqw(i,k))
1452#ifdef ISO
1453        do ixt=1,ntraciso
1454          omgbdxt(ixt,i,k)  = omgb(i,k)*(deltaxtw(ixt,i,k-1) - deltaxtw(ixt,i,k))
1455        enddo !do ixt=1,ntraciso
1456#ifdef ISOVERIF
1457      if (iso_eau.gt.0) then
1458        call iso_verif_egalite(deltaqw(i,k-1),deltaxtw(iso_eau,i,k-1),'wake 1460a')
1459        call iso_verif_egalite(deltaqw(i,k),deltaxtw(iso_eau,i,k),'wake 1460b')
1460        call iso_verif_egalite(omgbdq(i,k),omgbdxt(iso_eau,i,k),'wake 1460c')
1461      endif
1462#endif
1463#endif
1464        END IF
1465      END DO
1466    END DO
1467
1468    IF (prt_level>=10) THEN
1469      PRINT *, 'wake-4.3, th1(igout,k) ', (k,th1(igout,k), k=1,klev)
1470      PRINT *, 'wake-4.3, th2(igout,k) ', (k,th2(igout,k), k=1,klev)
1471      PRINT *, 'wake-4.3, dth(igout,k) ', (k,dth(igout,k), k=1,klev)
1472      PRINT *, 'wake-4.3, omgbdth(igout,k) ', (k,omgbdth(igout,k), k=1,klev)
1473    ENDIF
1474
1475    ! -----------------------------------------------------------------
1476    DO k = 1, klev-1
1477      DO i = 1, klon
1478        IF (wk_adv(i) .AND. k<=kupper(i)-1) THEN
1479          ! -----------------------------------------------------------------
1480
1481          ! Compute redistribution (advective) term
1482
1483          d_deltatw(i, k) = dtimesub/(ph(i,k)-ph(i,k+1))* &
1484            (rrd1*omg(i,k)*sigmaw(i)*d_th1(i,k) - &
1485             rrd2*omg(i,k+1)*(1.-sigmaw(i))*d_th2(i,k+1)- &
1486             (1.-alpha_up(i,k))*omgbdth(i,k)- &
1487             alpha_up(i,k+1)*omgbdth(i,k+1))*ppi(i, k)
1488!           print*,'d_deltatw=', k, d_deltatw(i,k)
1489
1490          d_deltaqw(i, k) = dtimesub/(ph(i,k)-ph(i,k+1))* &
1491            (rrd1*omg(i,k)*sigmaw(i)*d_q1(i,k)- &
1492             rrd2*omg(i,k+1)*(1.-sigmaw(i))*d_q2(i,k+1)- &
1493             (1.-alpha_up(i,k))*omgbdq(i,k)- &
1494             alpha_up(i,k+1)*omgbdq(i,k+1))
1495!           print*,'d_deltaqw=', k, d_deltaqw(i,k)
1496#ifdef ISO
1497        do ixt=1,ntraciso
1498          d_deltaxtw(ixt,i,k) = dtimesub/(Ph(i,k)-Ph(i,k+1))* &
1499            (rrd1*omg(i,k  )*sigmaw(i)     *d_xt1(ixt,i,k)- &
1500            rrd2*omg(i,k+1)*(1.-sigmaw(i))*d_xt2(ixt,i,k+1)- &
1501            (1.-alpha_up(i,k))*omgbdxt(ixt,i,k)- &
1502            alpha_up(i,k+1)*omgbdxt(ixt,i,k+1))
1503        enddo !do ixt=1,ntraciso
1504#ifdef ISOVERIF
1505      if (iso_eau.gt.0) then
1506        call iso_verif_egalite(d_q1(i,k),d_xt1(iso_eau,i,k),'wake 1502a')
1507        call iso_verif_egalite(d_q2(i,k),d_xt2(iso_eau,i,k),'wake 1502b')
1508        call iso_verif_egalite(omgbdq(i,k),omgbdxt(iso_eau,i,k),'wake 1502c')
1509        call iso_verif_egalite(omgbdq(i,k+1),omgbdxt(iso_eau,i,k+1),'wake 1502d')
1510        call iso_verif_egalite(d_deltaqw(i,k),d_deltaxtw(iso_eau,i,k),'wake 1502e')
1511      endif
1512#endif
1513#endif
1514
1515          ! and increment large scale tendencies
1516
1517
1518
1519
1520          ! C
1521          ! -----------------------------------------------------------------
1522          d_te(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_th1(i,k)- &
1523                                  rre2(i)*omg(i,k+1)*(1.-sigmaw(i))*d_th2(i,k+1))/ &
1524                                 (ph(i,k)-ph(i,k+1)) &
1525                                 -sigmaw(i)*(1.-sigmaw(i))*dth(i,k)*(omg(i,k)-omg(i,k+1))/ &
1526                                 (ph(i,k)-ph(i,k+1)) )*ppi(i, k)
1527
1528          d_qe(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_q1(i,k)- &
1529                                  rre2(i)*omg(i,k+1)*(1.-sigmaw(i))*d_q2(i,k+1))/ &
1530                                 (ph(i,k)-ph(i,k+1)) &
1531                                 -sigmaw(i)*(1.-sigmaw(i))*deltaqw(i,k)*(omg(i,k)-omg(i,k+1))/ &
1532                                 (ph(i,k)-ph(i,k+1)) )
1533#ifdef ISO
1534        do ixt=1,ntraciso
1535         d_xte(ixt,i,k) =  dtimesub*( &
1536             ( rre1(i)*omg(i,k  )*sigmaw(i)     *d_xt1(ixt,i,k) &
1537              -rre2(i)*omg(i,k+1)*(1.-sigmaw(i))*d_xt2(ixt,i,k+1) ) &
1538                    /(Ph(i,k)-Ph(i,k+1)) &
1539              -sigmaw(i)*(1.-sigmaw(i))*deltaxtw(ixt,i,k) &
1540              *(omg(i,k)-omg(i,k+1))/(Ph(i,k)-Ph(i,k+1)) &
1541              )
1542        enddo !do ixt=1,ntraciso
1543#endif
1544        ELSE IF (wk_adv(i) .AND. k==kupper(i)) THEN
1545          d_te(i, k) = dtimesub*(rre1(i)*omg(i,k)*sigmaw(i)*d_th1(i,k)/(ph(i,k)-ph(i,k+1)))*ppi(i, k)
1546
1547          d_qe(i, k) = dtimesub*(rre1(i)*omg(i,k)*sigmaw(i)*d_q1(i,k)/(ph(i,k)-ph(i,k+1)))
1548
1549#ifdef ISO
1550        do ixt=1,ntraciso
1551         d_xte(ixt,i,k) =  dtimesub*( &
1552            ( rre1(i)*omg(i,k  )*sigmaw(i)  *d_xt1(ixt,i,k)     &   
1553             /(Ph(i,k)-Ph(i,k+1))) &
1554                            )
1555        enddo !do ixt=1,ntraciso
1556#endif
1557        END IF
1558        ! cc
1559      END DO
1560    END DO
1561    ! ------------------------------------------------------------------
1562
1563    IF (prt_level>=10) THEN
1564      PRINT *, 'wake-4.3, d_deltatw(igout,k) ', (k,d_deltatw(igout,k), k=1,klev)
1565      PRINT *, 'wake-4.3, d_deltaqw(igout,k) ', (k,d_deltaqw(igout,k), k=1,klev)
1566    ENDIF
1567
1568    ! Increment state variables
1569!jyg<
1570    IF (iflag_wk_pop_dyn >= 1) THEN
1571      DO k = 1, klev
1572        DO i = 1, klon
1573          IF (wk_adv(i) .AND. k<=kupper(i)) THEN
1574            detr(i,k) = - d_sig_death(i) - d_sig_col(i)     
1575            entr(i,k) = d_sig_gen(i)
1576          ENDIF
1577        ENDDO
1578      ENDDO
1579      ELSE  ! (iflag_wk_pop_dyn >= 1)
1580      DO k = 1, klev
1581        DO i = 1, klon
1582          IF (wk_adv(i) .AND. k<=kupper(i)) THEN
1583            detr(i, k) = 0.
1584   
1585            entr(i, k) = 0.
1586          ENDIF
1587        ENDDO
1588      ENDDO
1589    ENDIF  ! (iflag_wk_pop_dyn >= 1)
1590
1591   
1592
1593    DO k = 1, klev
1594      DO i = 1, klon
1595        ! cc nrlmd       IF( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN
1596        IF (wk_adv(i) .AND. k<=kupper(i)) THEN
1597          ! cc
1598
1599
1600
1601          ! Coefficient de répartition
1602
1603          crep(i, k) = crep_sol*(ph(i,kupper(i))-ph(i,k))/ &
1604            (ph(i,kupper(i))-ph(i,1))
1605          crep(i, k) = crep(i, k) + crep_upper*(ph(i,1)-ph(i,k))/ &
1606            (p(i,1)-ph(i,kupper(i)))
1607
1608
1609          ! Reintroduce compensating subsidence term.
1610
1611          ! dtKE(k)=(dtdwn(k)*Crep(k))/sigmaw
1612          ! dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k))
1613          ! .                   /(1-sigmaw)
1614          ! dqKE(k)=(dqdwn(k)*Crep(k))/sigmaw
1615          ! dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k))
1616          ! .                   /(1-sigmaw)
1617
1618          ! dtKE(k)=(dtdwn(k)*Crep(k)+(1-Crep(k))*dta(k))/sigmaw
1619          ! dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k)*Crep(k))
1620          ! .                   /(1-sigmaw)
1621          ! dqKE(k)=(dqdwn(k)*Crep(k)+(1-Crep(k))*dqa(k))/sigmaw
1622          ! dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k)*Crep(k))
1623          ! .                   /(1-sigmaw)
1624
1625          dtke(i, k) = (dtdwn(i,k)/sigmaw(i)-dta(i,k)/(1.-sigmaw(i)))
1626          dqke(i, k) = (dqdwn(i,k)/sigmaw(i)-dqa(i,k)/(1.-sigmaw(i)))
1627#ifdef ISO
1628        do ixt=1,ntraciso
1629          dxtke(ixt,i,k)=(dxtdwn(ixt,i,k)/sigmaw(i) - dxta(ixt,i,k) &
1630               /(1.-sigmaw(i)))
1631        enddo !do ixt=1,ntraciso
1632#ifdef ISOVERIF
1633      if (iso_eau.gt.0) then
1634        call iso_verif_egalite(dqke(i,k),dxtKE(iso_eau,i,k),'wake 1621a')
1635        call iso_verif_egalite(dqdwn(i,k),dxtdwn(iso_eau,i,k),'wake 1621b')
1636        call iso_verif_egalite(dqa(i,k),dxta(iso_eau,i,k),'wake 1621c')
1637        call iso_verif_egalite(d_deltaqw(i,k),d_deltaxtw(iso_eau,i,k),'wake 1621d')
1638      endif
1639#endif
1640#endif
1641          ! print*,'dtKE= ',dtKE(i,k),' dqKE= ',dqKE(i,k)
1642
1643!
1644
1645          ! cc nrlmd          Prise en compte du taux de mortalité
1646          ! cc               Définitions de entr, detr
1647!jyg<
1648!!            detr(i, k) = 0.
1649!!   
1650!!            entr(i, k) = detr(i, k) + gfl(i)*cstar(i) + &
1651!!              sigmaw(i)*(1.-sigmaw(i))*dp_deltomg(i, k)
1652!!
1653            entr(i, k) = entr(i,k) + gfl(i)*cstar(i) + &
1654                         sigmaw(i)*(1.-sigmaw(i))*dp_deltomg(i, k)   
1655!>jyg
1656            spread(i, k) = (entr(i,k)-detr(i,k))/sigmaw(i)
1657
1658          ! cc        spread(i,k) =
1659          ! (1.-sigmaw(i))*dp_deltomg(i,k)+gfl(i)*Cstar(i)/
1660          ! cc     $  sigmaw(i)
1661
1662
1663          ! ajout d'un effet onde de gravité -Tgw(k)*deltatw(k) 03/02/06 YU
1664          ! Jingmei
1665
1666          ! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltat_gw(i,k),
1667          ! &  Tgw(i,k),deltatw(i,k)
1668          d_deltat_gw(i, k) = d_deltat_gw(i, k) - tgw(i, k)*deltatw(i, k)* &
1669            dtimesub
1670          ! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltatw(i,k)
1671          ff(i) = d_deltatw(i, k)/dtimesub
1672
1673          ! Sans GW
1674
1675          ! deltatw(k)=deltatw(k)+dtimesub*(ff+dtKE(k)-spread(k)*deltatw(k))
1676
1677          ! GW formule 1
1678
1679          ! deltatw(k) = deltatw(k)+dtimesub*
1680          ! $         (ff+dtKE(k) - spread(k)*deltatw(k)-Tgw(k)*deltatw(k))
1681
1682          ! GW formule 2
1683
1684          IF (dtimesub*tgw(i,k)<1.E-10) THEN
1685            d_deltatw(i, k) = dtimesub*(ff(i)+dtke(i,k) - &
1686               entr(i,k)*deltatw(i,k)/sigmaw(i) - &
1687               (death_rate(i)*sigmaw(i)+detr(i,k))*deltatw(i,k)/(1.-sigmaw(i)) - & ! cc
1688               tgw(i,k)*deltatw(i,k) )
1689          ELSE
1690            d_deltatw(i, k) = 1/tgw(i, k)*(1-exp(-dtimesub*tgw(i,k)))* &
1691               (ff(i)+dtke(i,k) - &
1692                entr(i,k)*deltatw(i,k)/sigmaw(i) - &
1693                (death_rate(i)*sigmaw(i)+detr(i,k))*deltatw(i,k)/(1.-sigmaw(i)) - &
1694                tgw(i,k)*deltatw(i,k) )
1695          END IF
1696
1697          dth(i, k) = deltatw(i, k)/ppi(i, k)
1698
1699          gg(i) = d_deltaqw(i, k)/dtimesub
1700
1701          d_deltaqw(i, k) = dtimesub*(gg(i)+dqke(i,k) - &
1702            entr(i,k)*deltaqw(i,k)/sigmaw(i) - &
1703            (death_rate(i)*sigmaw(i)+detr(i,k))*deltaqw(i,k)/(1.-sigmaw(i)))
1704#ifdef ISO
1705        do ixt=1,ntraciso
1706          gg(i)=d_deltaxtw(ixt,i,k)/dtimesub
1707          d_deltaxtw(ixt,i,k) = dtimesub*(gg(i) + dxtKE(ixt,i,k) - &
1708             entr(i,k)*deltaxtw(ixt,i,k)/sigmaw(i) - &
1709             (death_rate(i)*sigmaw(i)+detr(i,k))*deltaxtw(ixt,i,k)/(1.-sigmaw(i)))
1710        enddo !do ixt=1,ntraciso
1711#ifdef ISOVERIF
1712      if (iso_eau.gt.0) then
1713        call iso_verif_egalite(dqke(i,k),dxtKE(iso_eau,i,k),'wake 1692a')
1714        call iso_verif_egalite(deltaqw(i,k),deltaxtw(iso_eau,i,k),'wake 1692b')
1715        call iso_verif_egalite(d_deltaqw(i,k),d_deltaxtw(iso_eau,i,k),'wake 1692c')
1716      endif
1717#endif
1718#endif
1719          ! cc
1720
1721          ! cc nrlmd
1722          ! cc       d_deltatw2(i,k)=d_deltatw2(i,k)+d_deltatw(i,k)
1723          ! cc       d_deltaqw2(i,k)=d_deltaqw2(i,k)+d_deltaqw(i,k)
1724          ! cc
1725        END IF
1726      END DO
1727    END DO
1728
1729#ifdef ISO
1730#ifdef ISOVERIF
1731      if (iso_eau.gt.0) then
1732             call iso_verif_egalite_vect2D(d_deltaxtw,d_deltaqw, &
1733                 'wake 1359',ntraciso,klon,klev)
1734      endif     
1735#endif         
1736#endif 
1737
1738    ! Scale tendencies so that water vapour remains positive in w and x.
1739
1740    CALL wake_vec_modulation(klon, klev, wk_adv, epsilon, qe, d_qe, deltaqw, &
1741      d_deltaqw, sigmaw, d_sigmaw, alpha)
1742
1743    ! cc nrlmd
1744    ! c      print*,'alpha'
1745    ! c      do i=1,klon
1746    ! c         print*,alpha(i)
1747    ! c      end do
1748    ! cc
1749    DO k = 1, klev
1750      DO i = 1, klon
1751        IF (wk_adv(i) .AND. k<=kupper(i)) THEN
1752          d_te(i, k) = alpha(i)*d_te(i, k)
1753          d_qe(i, k) = alpha(i)*d_qe(i, k)
1754          d_deltatw(i, k) = alpha(i)*d_deltatw(i, k)
1755          d_deltaqw(i, k) = alpha(i)*d_deltaqw(i, k)
1756          d_deltat_gw(i, k) = alpha(i)*d_deltat_gw(i, k)
1757#ifdef ISO
1758        do ixt=1,ntraciso
1759          d_xte(ixt,i,k)=alpha(i)*d_xte(ixt,i,k)
1760          d_deltaxtw(ixt,i,k)=alpha(i)*d_deltaxtw(ixt,i,k)
1761        enddo !do ixt=1,ntraciso
1762#endif
1763        END IF
1764      END DO
1765    END DO
1766    DO i = 1, klon
1767      IF (wk_adv(i)) THEN
1768        d_sigmaw(i) = alpha(i)*d_sigmaw(i)
1769      END IF
1770    END DO
1771
1772    ! Update large scale variables and wake variables
1773    ! IM 060208 manque DO i + remplace DO k=1,kupper(i)
1774    ! IM 060208     DO k = 1,kupper(i)
1775    DO k = 1, klev
1776      DO i = 1, klon
1777        IF (wk_adv(i) .AND. k<=kupper(i)) THEN
1778          dtls(i, k) = dtls(i, k) + d_te(i, k)
1779          dqls(i, k) = dqls(i, k) + d_qe(i, k)
1780#ifdef ISO
1781        do ixt=1,ntraciso
1782          dxtls(ixt,i,k)=dxtls(ixt,i,k)+d_xte(ixt,i,k)
1783        enddo !do ixt=1,ntraciso
1784#endif
1785          ! cc nrlmd
1786          d_deltatw2(i, k) = d_deltatw2(i, k) + d_deltatw(i, k)
1787          d_deltaqw2(i, k) = d_deltaqw2(i, k) + d_deltaqw(i, k)
1788#ifdef ISO
1789        do ixt=1,ntraciso
1790          d_deltaxtw2(ixt,i,k)=d_deltaxtw2(ixt,i,k)+d_deltaxtw(ixt,i,k)
1791        enddo !do ixt=1,ntraciso
1792#endif
1793          ! cc
1794        END IF
1795      END DO
1796    END DO
1797
1798
1799#ifdef ISO
1800#ifdef ISOVERIF
1801      if (iso_eau.gt.0) then
1802        DO k= 1,klev
1803          DO i = 1,klon
1804              call iso_verif_egalite_choix(dxtls(iso_eau,i,k), &
1805                dqls(i,k),'wake 1379',errmax,errmaxrel)
1806          enddo ! DO i = 1,klon     
1807        enddo ! DO k= 1,klev
1808      endif !if (iso_eau.gt.0) then
1809#endif
1810#endif
1811
1812
1813    DO k = 1, klev
1814      DO i = 1, klon
1815        IF (wk_adv(i) .AND. k<=kupper(i)) THEN
1816          te(i, k) = te0(i, k) + dtls(i, k)
1817          qe(i, k) = qe0(i, k) + dqls(i, k)
1818          the(i, k) = te(i, k)/ppi(i, k)
1819          deltatw(i, k) = deltatw(i, k) + d_deltatw(i, k)
1820          deltaqw(i, k) = deltaqw(i, k) + d_deltaqw(i, k)
1821          dth(i, k) = deltatw(i, k)/ppi(i, k)
1822          ! c      print*,'k,qx,qw',k,qe(i,k)-sigmaw(i)*deltaqw(i,k)
1823          ! c     $        ,qe(i,k)+(1-sigmaw(i))*deltaqw(i,k)
1824#ifdef ISO
1825        do ixt=1,ntraciso
1826          xte(ixt,i,k) = xte0(ixt,i,k) + dxtls(ixt,i,k)
1827          deltaxtw(ixt,i,k) = deltaxtw(ixt,i,k)+d_deltaxtw(ixt,i,k)
1828        enddo !do ixt=1,ntraciso
1829#endif
1830        END IF
1831      END DO
1832    END DO
1833!
1834    DO i = 1, klon
1835      IF (wk_adv(i)) THEN
1836        sigmaw(i) = sigmaw(i) + d_sigmaw(i)
1837        d_sigmaw2(i) = d_sigmaw2(i) + d_sigmaw(i)
1838      END IF
1839    END DO
1840
1841#ifdef ISO
1842#ifdef ISOVERIF
1843      write(*,*) 'wake 1859'
1844      if (iso_eau.gt.0) then
1845        DO k= 1,klev
1846          DO i = 1,klon
1847              call iso_verif_egalite_choix(xte(iso_eau,i,k), &
1848                qe(i,k),'wake 1379',errmax,errmaxrel)
1849          enddo ! DO i = 1,klon     
1850        enddo ! DO k= 1,klev
1851      endif !if (iso_eau.gt.0) then     
1852      if (iso_hdo.gt.0) then
1853      call iso_verif_aberrant_enc_vect2D( &
1854            xte,qe, &
1855            'wake 1456, xte apres modifs',ntraciso,klon,klev)
1856!      call iso_verif_aberrant_enc_vect2D_ns(
1857!     :       deltaxtw,deltaqw,
1858!     :       'wake 1518, deltaqw apres modifs',ntraciso,klon,klev)
1859      endif
1860#endif
1861#endif
1862
1863!jyg<
1864    IF (iflag_wk_pop_dyn >= 1) THEN
1865      DO i = 1, klon
1866        IF (wk_adv(i)) THEN
1867          awdens(i) = awdens(i) + d_awdens(i)
1868          wdens(i) = wdens(i) + d_wdens(i)
1869          d_awdens2(i) = d_awdens2(i) + d_awdens(i)
1870          d_wdens2(i) = d_wdens2(i) + d_wdens(i)
1871        END IF
1872      END DO
1873      DO i = 1, klon
1874        IF (wk_adv(i)) THEN
1875          wdens_targ = max(wdens(i),wdensmin)
1876          d_wdens2(i) = d_wdens2(i) + wdens_targ - wdens(i)
1877          wdens(i) = wdens_targ
1878!
1879          wdens_targ = min( max(awdens(i),0.), wdens(i) )
1880          d_awdens2(i) = d_awdens2(i) + wdens_targ - awdens(i)
1881          awdens(i) = wdens_targ
1882        END IF
1883      END DO
1884      DO i = 1, klon
1885        IF (wk_adv(i)) THEN
1886          sigmaw_targ = max(sigmaw(i),sigmad)
1887          d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
1888          sigmaw(i) = sigmaw_targ
1889        END IF
1890      END DO
1891    ENDIF  ! (iflag_wk_pop_dyn >= 1)
1892!>jyg
1893
1894
1895    ! Determine Ptop from buoyancy integral
1896    ! ---------------------------------------
1897
1898    ! -     1/ Pressure of the level where dth changes sign.
1899
1900    DO i = 1, klon
1901      IF (wk_adv(i)) THEN
1902        ptop_provis(i) = ph(i, 1)
1903      END IF
1904    END DO
1905
1906    DO k = 2, klev
1907      DO i = 1, klon
1908        IF (wk_adv(i) .AND. ptop_provis(i)==ph(i,1) .AND. &
1909            dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
1910          ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) - &
1911                            (dth(i,k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
1912        END IF
1913      END DO
1914    END DO
1915
1916    ! -     2/ dth integral
1917
1918    DO i = 1, klon
1919      IF (wk_adv(i)) THEN !!! nrlmd
1920        sum_dth(i) = 0.
1921        dthmin(i) = -delta_t_min
1922        z(i) = 0.
1923      END IF
1924    END DO
1925
1926    DO k = 1, klev
1927      DO i = 1, klon
1928        IF (wk_adv(i)) THEN
1929          dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*rg)
1930          IF (dz(i)>0) THEN
1931            z(i) = z(i) + dz(i)
1932            sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
1933            dthmin(i) = amin1(dthmin(i), dth(i,k))
1934          END IF
1935        END IF
1936      END DO
1937    END DO
1938
1939    ! -     3/ height of triangle with area= sum_dth and base = dthmin
1940
1941    DO i = 1, klon
1942      IF (wk_adv(i)) THEN
1943        hw(i) = 2.*sum_dth(i)/amin1(dthmin(i), -0.5)
1944        hw(i) = amax1(hwmin, hw(i))
1945      END IF
1946    END DO
1947
1948    ! -     4/ now, get Ptop
1949
1950    DO i = 1, klon
1951      IF (wk_adv(i)) THEN !!! nrlmd
1952        ktop(i) = 0
1953        z(i) = 0.
1954      END IF
1955    END DO
1956
1957    DO k = 1, klev
1958      DO i = 1, klon
1959        IF (wk_adv(i)) THEN
1960          dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg), hw(i)-z(i))
1961          IF (dz(i)>0) THEN
1962            z(i) = z(i) + dz(i)
1963            ptop(i) = ph(i, k) - rho(i, k)*rg*dz(i)
1964            ktop(i) = k
1965          END IF
1966        END IF
1967      END DO
1968    END DO
1969
1970    ! 4.5/Correct ktop and ptop
1971
1972    DO i = 1, klon
1973      IF (wk_adv(i)) THEN
1974        ptop_new(i) = ptop(i)
1975      END IF
1976    END DO
1977
1978    DO k = klev, 2, -1
1979      DO i = 1, klon
1980        ! IM v3JYG; IF (k .GE. ktop(i)
1981        IF (wk_adv(i) .AND. k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. &
1982            dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
1983          ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) - &
1984                         (dth(i,k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
1985        END IF
1986      END DO
1987    END DO
1988
1989
1990    DO i = 1, klon
1991      IF (wk_adv(i)) THEN
1992        ptop(i) = ptop_new(i)
1993      END IF
1994    END DO
1995
1996    DO k = klev, 1, -1
1997      DO i = 1, klon
1998        IF (wk_adv(i)) THEN !!! nrlmd
1999          IF (ph(i,k+1)<ptop(i)) ktop(i) = k
2000        END IF
2001      END DO
2002    END DO
2003
2004    ! 5/ Set deltatw & deltaqw to 0 above kupper
2005
2006    DO k = 1, klev
2007      DO i = 1, klon
2008        IF (wk_adv(i) .AND. k>=kupper(i)) THEN
2009          deltatw(i, k) = 0.
2010          deltaqw(i, k) = 0.
2011          d_deltatw2(i,k) = -deltatw0(i,k)
2012          d_deltaqw2(i,k) = -deltaqw0(i,k)
2013#ifdef ISO
2014        do ixt=1,ntraciso
2015          deltaxtw(ixt,i,k) = 0.
2016          d_deltaxtw2(ixt,i,k) = -deltaxtw0(ixt,i,k)
2017        enddo !do ixt=1,ntraciso
2018#endif
2019        END IF
2020      END DO
2021    END DO
2022
2023
2024    ! -------------Cstar computation---------------------------------
2025    DO i = 1, klon
2026      IF (wk_adv(i)) THEN !!! nrlmd
2027        sum_thu(i) = 0.
2028        sum_tu(i) = 0.
2029        sum_qu(i) = 0.
2030        sum_thvu(i) = 0.
2031        sum_dth(i) = 0.
2032        sum_dq(i) = 0.
2033        sum_rho(i) = 0.
2034        sum_dtdwn(i) = 0.
2035        sum_dqdwn(i) = 0.
2036
2037        av_thu(i) = 0.
2038        av_tu(i) = 0.
2039        av_qu(i) = 0.
2040        av_thvu(i) = 0.
2041        av_dth(i) = 0.
2042        av_dq(i) = 0.
2043        av_rho(i) = 0.
2044        av_dtdwn(i) = 0.
2045        av_dqdwn(i) = 0.
2046      END IF
2047    END DO
2048
2049    ! Integrals (and wake top level number)
2050    ! --------------------------------------
2051
2052    ! Initialize sum_thvu to 1st level virt. pot. temp.
2053
2054    DO i = 1, klon
2055      IF (wk_adv(i)) THEN !!! nrlmd
2056        z(i) = 1.
2057        dz(i) = 1.
2058        sum_thvu(i) = thu(i, 1)*(1.+epsim1*qu(i,1))*dz(i)
2059        sum_dth(i) = 0.
2060      END IF
2061    END DO
2062
2063    DO k = 1, klev
2064      DO i = 1, klon
2065        IF (wk_adv(i)) THEN !!! nrlmd
2066          dz(i) = -(max(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
2067          IF (dz(i)>0) THEN
2068            z(i) = z(i) + dz(i)
2069            sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i)
2070            sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i)
2071            sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i)
2072            sum_thvu(i) = sum_thvu(i) + thu(i, k)*(1.+epsim1*qu(i,k))*dz(i)
2073            sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
2074            sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
2075            sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i)
2076            sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
2077            sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
2078          END IF
2079        END IF
2080      END DO
2081    END DO
2082
2083    DO i = 1, klon
2084      IF (wk_adv(i)) THEN !!! nrlmd
2085        hw0(i) = z(i)
2086      END IF
2087    END DO
2088
2089
2090    ! - WAPE and mean forcing computation
2091    ! ---------------------------------------
2092
2093    ! ---------------------------------------
2094
2095    ! Means
2096
2097    DO i = 1, klon
2098      IF (wk_adv(i)) THEN !!! nrlmd
2099        av_thu(i) = sum_thu(i)/hw0(i)
2100        av_tu(i) = sum_tu(i)/hw0(i)
2101        av_qu(i) = sum_qu(i)/hw0(i)
2102        av_thvu(i) = sum_thvu(i)/hw0(i)
2103        av_dth(i) = sum_dth(i)/hw0(i)
2104        av_dq(i) = sum_dq(i)/hw0(i)
2105        av_rho(i) = sum_rho(i)/hw0(i)
2106        av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
2107        av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
2108
2109        wape(i) = -rg*hw0(i)*(av_dth(i)+epsim1*(av_thu(i)*av_dq(i) + &
2110                              av_dth(i)*av_qu(i)+av_dth(i)*av_dq(i)))/av_thvu(i)
2111      END IF
2112    END DO
2113
2114    ! Filter out bad wakes
2115
2116    DO k = 1, klev
2117      DO i = 1, klon
2118        IF (wk_adv(i)) THEN !!! nrlmd
2119          IF (wape(i)<0.) THEN
2120            deltatw(i, k) = 0.
2121            deltaqw(i, k) = 0.
2122            dth(i, k) = 0.
2123            d_deltatw2(i,k) = -deltatw0(i,k)
2124            d_deltaqw2(i,k) = -deltaqw0(i,k)
2125#ifdef ISO
2126        do ixt=1,ntraciso
2127          deltaxtw(ixt,i,k) = 0.
2128          d_deltaxtw2(ixt,i,k) = -deltaxtw0(ixt,i,k)
2129        enddo !do ixt=1,ntraciso
2130#endif
2131          END IF
2132        END IF
2133      END DO
2134    END DO
2135
2136    DO i = 1, klon
2137      IF (wk_adv(i)) THEN !!! nrlmd
2138        IF (wape(i)<0.) THEN
2139          wape(i) = 0.
2140          cstar(i) = 0.
2141          hw(i) = hwmin
2142!jyg<
2143!!          sigmaw(i) = max(sigmad, sigd_con(i))
2144          sigmaw_targ = max(sigmad, sigd_con(i))
2145          d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
2146          sigmaw(i) = sigmaw_targ
2147!>jyg
2148          fip(i) = 0.
2149          gwake(i) = .FALSE.
2150        ELSE
2151          cstar(i) = stark*sqrt(2.*wape(i))
2152          gwake(i) = .TRUE.
2153        END IF
2154      END IF
2155    END DO
2156
2157  END DO ! end sub-timestep loop
2158
2159  IF (prt_level>=10) THEN
2160    PRINT *, 'wake-5, sigmaw(igout), cstar(igout), wape(igout), ptop(igout) ', &
2161                      sigmaw(igout), cstar(igout), wape(igout), ptop(igout)
2162  ENDIF
2163
2164
2165  ! ----------------------------------------------------------
2166  ! Determine wake final state; recompute wape, cstar, ktop;
2167  ! filter out bad wakes.
2168  ! ----------------------------------------------------------
2169
2170  ! 2.1 - Undisturbed area and Wake integrals
2171  ! ---------------------------------------------------------
2172
2173  DO i = 1, klon
2174    ! cc nrlmd       if (wk_adv(i)) then !!! nrlmd
2175    IF (ok_qx_qw(i)) THEN
2176      ! cc
2177      z(i) = 0.
2178      sum_thu(i) = 0.
2179      sum_tu(i) = 0.
2180      sum_qu(i) = 0.
2181      sum_thvu(i) = 0.
2182      sum_dth(i) = 0.
2183      sum_half_dth(i) = 0.
2184      sum_dq(i) = 0.
2185      sum_rho(i) = 0.
2186      sum_dtdwn(i) = 0.
2187      sum_dqdwn(i) = 0.
2188
2189      av_thu(i) = 0.
2190      av_tu(i) = 0.
2191      av_qu(i) = 0.
2192      av_thvu(i) = 0.
2193      av_dth(i) = 0.
2194      av_dq(i) = 0.
2195      av_rho(i) = 0.
2196      av_dtdwn(i) = 0.
2197      av_dqdwn(i) = 0.
2198
2199      dthmin(i) = -delta_t_min
2200    END IF
2201  END DO
2202  ! Potential temperatures and humidity
2203  ! ----------------------------------------------------------
2204
2205  DO k = 1, klev
2206    DO i = 1, klon
2207      ! cc nrlmd       IF ( wk_adv(i)) THEN
2208      IF (ok_qx_qw(i)) THEN
2209        ! cc
2210        rho(i, k) = p(i, k)/(rd*te(i,k))
2211        IF (k==1) THEN
2212          rhoh(i, k) = ph(i, k)/(rd*te(i,k))
2213          zhh(i, k) = 0
2214        ELSE
2215          rhoh(i, k) = ph(i, k)*2./(rd*(te(i,k)+te(i,k-1)))
2216          zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*rg) + zhh(i, k-1)
2217        END IF
2218        the(i, k) = te(i, k)/ppi(i, k)
2219        thu(i, k) = (te(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k)
2220        tu(i, k) = te(i, k) - deltatw(i, k)*sigmaw(i)
2221        qu(i, k) = qe(i, k) - deltaqw(i, k)*sigmaw(i)
2222        rhow(i, k) = p(i, k)/(rd*(te(i,k)+deltatw(i,k)))
2223        dth(i, k) = deltatw(i, k)/ppi(i, k)
2224#ifdef ISO
2225        do ixt=1,ntraciso
2226          xtu(ixt,i,k)  =  xte(ixt,i,k) - deltaxtw(ixt,i,k)*sigmaw(i)
2227        enddo !do ixt=1,ntraciso
2228#endif
2229      END IF
2230    END DO
2231  END DO
2232
2233#ifdef ISO
2234#ifdef ISOVERIF
2235      if (iso_hdo.gt.0) then
2236      call iso_verif_aberrant_enc_vect2D( &
2237              xtu,qu, &
2238              'wake 1834, apres modifs',ntraciso,klon,klev)
2239      endif   
2240#endif
2241#endif   
2242
2243  ! Integrals (and wake top level number)
2244  ! -----------------------------------------------------------
2245
2246  ! Initialize sum_thvu to 1st level virt. pot. temp.
2247
2248  DO i = 1, klon
2249    ! cc nrlmd       IF ( wk_adv(i)) THEN
2250    IF (ok_qx_qw(i)) THEN
2251      ! cc
2252      z(i) = 1.
2253      dz(i) = 1.
2254      dz_half(i) = 1.
2255      sum_thvu(i) = thu(i, 1)*(1.+epsim1*qu(i,1))*dz(i)
2256      sum_dth(i) = 0.
2257    END IF
2258  END DO
2259
2260  DO k = 1, klev
2261    DO i = 1, klon
2262      ! cc nrlmd       IF ( wk_adv(i)) THEN
2263      IF (ok_qx_qw(i)) THEN
2264        ! cc
2265        dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
2266        dz_half(i) = -(amax1(ph(i,k+1),0.5*(ptop(i)+ph(i,1)))-ph(i,k))/(rho(i,k)*rg)
2267        IF (dz(i)>0) THEN
2268          z(i) = z(i) + dz(i)
2269          sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i)
2270          sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i)
2271          sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i)
2272          sum_thvu(i) = sum_thvu(i) + thu(i, k)*(1.+epsim1*qu(i,k))*dz(i)
2273          sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
2274          sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
2275          sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i)
2276          sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
2277          sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
2278!
2279          dthmin(i) = min(dthmin(i), dth(i,k))
2280        END IF
2281        IF (dz_half(i)>0) THEN
2282          sum_half_dth(i) = sum_half_dth(i) + dth(i, k)*dz_half(i)
2283        END IF
2284      END IF
2285    END DO
2286  END DO
2287
2288  DO i = 1, klon
2289    ! cc nrlmd       IF ( wk_adv(i)) THEN
2290    IF (ok_qx_qw(i)) THEN
2291      ! cc
2292      hw0(i) = z(i)
2293    END IF
2294  END DO
2295
2296  ! - WAPE and mean forcing computation
2297  ! -------------------------------------------------------------
2298
2299  ! Means
2300
2301  DO i = 1, klon
2302    ! cc nrlmd       IF ( wk_adv(i)) THEN
2303    IF (ok_qx_qw(i)) THEN
2304      ! cc
2305      av_thu(i) = sum_thu(i)/hw0(i)
2306      av_tu(i) = sum_tu(i)/hw0(i)
2307      av_qu(i) = sum_qu(i)/hw0(i)
2308      av_thvu(i) = sum_thvu(i)/hw0(i)
2309      av_dth(i) = sum_dth(i)/hw0(i)
2310      av_dq(i) = sum_dq(i)/hw0(i)
2311      av_rho(i) = sum_rho(i)/hw0(i)
2312      av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
2313      av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
2314
2315      wape2(i) = -rg*hw0(i)*(av_dth(i)+epsim1*(av_thu(i)*av_dq(i) + &
2316                             av_dth(i)*av_qu(i)+av_dth(i)*av_dq(i)))/av_thvu(i)
2317    END IF
2318  END DO
2319
2320
2321
2322  ! Prognostic variable update
2323  ! ------------------------------------------------------------
2324
2325  ! Filter out bad wakes
2326
2327  IF (iflag_wk_check_trgl>=1) THEN
2328    ! Check triangular shape of dth profile
2329    DO i = 1, klon
2330      IF (ok_qx_qw(i)) THEN
2331        !! print *,'wake, hw0(i), dthmin(i) ', hw0(i), dthmin(i)
2332        !! print *,'wake, 2.*sum_dth(i)/(hw0(i)*dthmin(i)) ', &
2333        !!                2.*sum_dth(i)/(hw0(i)*dthmin(i))
2334        !! print *,'wake, sum_half_dth(i), sum_dth(i) ', &
2335        !!                sum_half_dth(i), sum_dth(i)
2336        IF ((hw0(i) < 1.) .or. (dthmin(i) >= -delta_t_min) ) THEN
2337          wape2(i) = -1.
2338          !! print *,'wake, rej 1'
2339        ELSE IF (iflag_wk_check_trgl==1.AND.abs(2.*sum_dth(i)/(hw0(i)*dthmin(i)) - 1.) > 0.5) THEN
2340          wape2(i) = -1.
2341          !! print *,'wake, rej 2'
2342        ELSE IF (abs(sum_half_dth(i)) < 0.5*abs(sum_dth(i)) ) THEN
2343          wape2(i) = -1.
2344          !! print *,'wake, rej 3'
2345        END IF
2346      END IF
2347    END DO
2348  END IF
2349
2350
2351  DO k = 1, klev
2352    DO i = 1, klon
2353      ! cc nrlmd        IF ( wk_adv(i) .AND. wape2(i) .LT. 0.) THEN
2354      IF (ok_qx_qw(i) .AND. wape2(i)<0.) THEN
2355        ! cc
2356        deltatw(i, k) = 0.
2357        deltaqw(i, k) = 0.
2358        dth(i, k) = 0.
2359        d_deltatw2(i,k) = -deltatw0(i,k)
2360        d_deltaqw2(i,k) = -deltaqw0(i,k)
2361#ifdef ISO
2362        do ixt=1,ntraciso
2363          deltaxtw(ixt,i,k) = 0.
2364          d_deltaxtw2(ixt,i,k) = -deltaxtw0(ixt,i,k)
2365        enddo !do ixt=1,ntraciso
2366#endif
2367      END IF
2368    END DO
2369  END DO
2370
2371
2372  DO i = 1, klon
2373    ! cc nrlmd       IF ( wk_adv(i)) THEN
2374    IF (ok_qx_qw(i)) THEN
2375      ! cc
2376      IF (wape2(i)<0.) THEN
2377        wape2(i) = 0.
2378        cstar2(i) = 0.
2379        hw(i) = hwmin
2380!jyg<
2381!!      sigmaw(i) = amax1(sigmad, sigd_con(i))
2382      sigmaw_targ = max(sigmad, sigd_con(i))
2383      d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
2384      sigmaw(i) = sigmaw_targ
2385!>jyg
2386        fip(i) = 0.
2387        gwake(i) = .FALSE.
2388      ELSE
2389        IF (prt_level>=10) PRINT *, 'wape2>0'
2390        cstar2(i) = stark*sqrt(2.*wape2(i))
2391        gwake(i) = .TRUE.
2392      END IF
2393    END IF
2394  END DO
2395
2396  DO i = 1, klon
2397    ! cc nrlmd       IF ( wk_adv(i)) THEN
2398    IF (ok_qx_qw(i)) THEN
2399      ! cc
2400      ktopw(i) = ktop(i)
2401    END IF
2402  END DO
2403
2404  DO i = 1, klon
2405    ! cc nrlmd       IF ( wk_adv(i)) THEN
2406    IF (ok_qx_qw(i)) THEN
2407      ! cc
2408      IF (ktopw(i)>0 .AND. gwake(i)) THEN
2409
2410        ! jyg1     Utilisation d'un h_efficace constant ( ~ feeding layer)
2411        ! cc       heff = 600.
2412        ! Utilisation de la hauteur hw
2413        ! c       heff = 0.7*hw
2414        heff(i) = hw(i)
2415
2416        fip(i) = 0.5*rho(i, ktopw(i))*cstar2(i)**3*heff(i)*2* &
2417          sqrt(sigmaw(i)*wdens(i)*3.14)
2418        fip(i) = alpk*fip(i)
2419        ! jyg2
2420      ELSE
2421        fip(i) = 0.
2422      END IF
2423    END IF
2424  END DO
2425
2426  ! Limitation de sigmaw
2427
2428  ! cc nrlmd
2429  ! DO i=1,klon
2430  ! IF (OK_qx_qw(i)) THEN
2431  ! IF (sigmaw(i).GE.sigmaw_max) sigmaw(i)=sigmaw_max
2432  ! ENDIF
2433  ! ENDDO
2434  ! cc
2435
2436  !jyg<
2437  IF (iflag_wk_pop_dyn >= 1) THEN
2438    DO i = 1, klon
2439      kill_wake(i) = ((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
2440          .NOT. ok_qx_qw(i) .OR. (wdens(i) < 2.*wdensmin)
2441    ENDDO
2442  ELSE  ! (iflag_wk_pop_dyn >= 1)
2443    DO i = 1, klon
2444      kill_wake(i) = ((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
2445          .NOT. ok_qx_qw(i)
2446    ENDDO
2447  ENDIF  ! (iflag_wk_pop_dyn >= 1)
2448  !>jyg
2449
2450  DO k = 1, klev
2451    DO i = 1, klon
2452!!jyg      IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
2453!!jyg          .NOT. ok_qx_qw(i)) THEN
2454      IF (kill_wake(i)) THEN
2455        ! cc
2456        dtls(i, k) = 0.
2457        dqls(i, k) = 0.
2458        deltatw(i, k) = 0.
2459        deltaqw(i, k) = 0.
2460        d_deltatw2(i,k) = -deltatw0(i,k)
2461        d_deltaqw2(i,k) = -deltaqw0(i,k)
2462#ifdef ISO
2463        do ixt=1,ntraciso
2464          dxtls(ixt,i,k) = 0.
2465          deltaxtw(ixt,i,k) = 0.
2466          d_deltaxtw2(ixt,i,k) = -deltaxtw0(ixt,i,k)
2467        enddo !do ixt=1,ntraciso
2468#endif
2469      END IF  ! (kill_wake(i))
2470    END DO
2471  END DO
2472
2473  DO i = 1, klon
2474!!jyg    IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
2475!!jyg        .NOT. ok_qx_qw(i)) THEN
2476      IF (kill_wake(i)) THEN
2477      ktopw(i) = 0
2478      wape(i) = 0.
2479      cstar(i) = 0.
2480!!jyg   Outside subroutine "Wake" hw, wdens and sigmaw are zero when there are no wakes
2481!!      hw(i) = hwmin                       !jyg
2482!!      sigmaw(i) = sigmad                  !jyg
2483      hw(i) = 0.                            !jyg
2484      fip(i) = 0.
2485!!      sigmaw(i) = 0.                        !jyg
2486      sigmaw_targ = 0.
2487      d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
2488      sigmaw(i) = sigmaw_targ
2489      IF (iflag_wk_pop_dyn >= 1) THEN
2490!!        awdens(i) = 0.
2491!!        wdens(i) = 0.
2492        wdens_targ = 0.
2493        d_wdens2(i) = wdens_targ - wdens(i)
2494        wdens(i) = wdens_targ
2495        wdens_targ = 0.
2496        d_awdens2(i) = wdens_targ - awdens(i)
2497        awdens(i) = wdens_targ
2498      ENDIF  ! (iflag_wk_pop_dyn >= 1)
2499    ELSE  ! (kill_wake(i))
2500      wape(i) = wape2(i)
2501      cstar(i) = cstar2(i)
2502    END IF  ! (kill_wake(i))
2503    ! c        print*,'wape wape2 ktopw OK_qx_qw =',
2504    ! c     $          wape(i),wape2(i),ktopw(i),OK_qx_qw(i)
2505  END DO
2506
2507  IF (prt_level>=10) THEN
2508    PRINT *, 'wake-6, wape wape2 ktopw OK_qx_qw =', &
2509                      wape(igout),wape2(igout),ktopw(igout),OK_qx_qw(igout)
2510  ENDIF
2511
2512
2513  ! -----------------------------------------------------------------
2514  ! Get back to tendencies per second
2515
2516  DO k = 1, klev
2517    DO i = 1, klon
2518
2519      ! cc nrlmd        IF ( wk_adv(i) .AND. k .LE. kupper(i)) THEN
2520!jyg<
2521!!      IF (ok_qx_qw(i) .AND. k<=kupper(i)) THEN
2522      IF (ok_qx_qw(i)) THEN
2523!>jyg
2524        ! cc
2525        dtls(i, k) = dtls(i, k)/dtime
2526        dqls(i, k) = dqls(i, k)/dtime
2527        d_deltatw2(i, k) = d_deltatw2(i, k)/dtime
2528        d_deltaqw2(i, k) = d_deltaqw2(i, k)/dtime
2529        d_deltat_gw(i, k) = d_deltat_gw(i, k)/dtime
2530        ! c      print*,'k,dqls,omg,entr,detr',k,dqls(i,k),omg(i,k),entr(i,k)
2531        ! c     $         ,death_rate(i)*sigmaw(i)
2532#ifdef ISO
2533        do ixt=1,ntraciso
2534          dxtls(ixt,i, k) = dxtls(ixt,i, k)/dtime
2535          d_deltaxtw2(ixt,i, k) = d_deltaxtw2(ixt,i, k)/dtime
2536        enddo
2537#endif
2538      END IF
2539    END DO
2540  END DO
2541!jyg<
2542  DO i = 1, klon
2543    d_sigmaw2(i) = d_sigmaw2(i)/dtime
2544    d_awdens2(i) = d_awdens2(i)/dtime
2545    d_wdens2(i) = d_wdens2(i)/dtime
2546  ENDDO
2547!>jyg
2548
2549
2550
2551  RETURN
2552END SUBROUTINE wake
2553
2554SUBROUTINE wake_vec_modulation(nlon, nl, wk_adv, epsilon, qe, d_qe, deltaqw, &
2555    d_deltaqw, sigmaw, d_sigmaw, alpha)
2556  ! ------------------------------------------------------
2557  ! D\'etermination du coefficient alpha tel que les tendances
2558  ! corriges alpha*d_G, pour toutes les grandeurs G, correspondent
2559  ! a une humidite positive dans la zone (x) et dans la zone (w).
2560  ! ------------------------------------------------------
2561  IMPLICIT NONE
2562
2563  ! Input
2564  REAL qe(nlon, nl), d_qe(nlon, nl)
2565  REAL deltaqw(nlon, nl), d_deltaqw(nlon, nl)
2566  REAL sigmaw(nlon), d_sigmaw(nlon)
2567  LOGICAL wk_adv(nlon)
2568  INTEGER nl, nlon
2569  ! Output
2570  REAL alpha(nlon)
2571  ! Internal variables
2572  REAL zeta(nlon, nl)
2573  REAL alpha1(nlon)
2574  REAL x, a, b, c, discrim
2575  REAL epsilon
2576  ! DATA epsilon/1.e-15/
2577  INTEGER i,k
2578
2579  DO k = 1, nl
2580    DO i = 1, nlon
2581      IF (wk_adv(i)) THEN
2582        IF ((deltaqw(i,k)+d_deltaqw(i,k))>=0.) THEN
2583          zeta(i, k) = 0.
2584        ELSE
2585          zeta(i, k) = 1.
2586        END IF
2587      END IF
2588    END DO
2589    DO i = 1, nlon
2590      IF (wk_adv(i)) THEN
2591        x = qe(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + d_qe(i, k) + &
2592          (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - d_sigmaw(i) * &
2593          (deltaqw(i,k)+d_deltaqw(i,k))
2594        a = -d_sigmaw(i)*d_deltaqw(i, k)
2595        b = d_qe(i, k) + (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - &
2596          deltaqw(i, k)*d_sigmaw(i)
2597        c = qe(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + epsilon
2598        discrim = b*b - 4.*a*c
2599        ! print*, 'x, a, b, c, discrim', x, a, b, c, discrim
2600        IF (a+b>=0.) THEN !! Condition suffisante pour la positivité de ovap
2601          alpha1(i) = 1.
2602        ELSE
2603          IF (x>=0.) THEN
2604            alpha1(i) = 1.
2605          ELSE
2606            IF (a>0.) THEN
2607              alpha1(i) = 0.9*min( (2.*c)/(-b+sqrt(discrim)),  &
2608                                   (-b+sqrt(discrim))/(2.*a) )
2609            ELSE IF (a==0.) THEN
2610              alpha1(i) = 0.9*(-c/b)
2611            ELSE
2612              ! print*,'a,b,c discrim',a,b,c discrim
2613              alpha1(i) = 0.9*max( (2.*c)/(-b+sqrt(discrim)),  &
2614                                   (-b+sqrt(discrim))/(2.*a))
2615            END IF
2616          END IF
2617        END IF
2618        alpha(i) = min(alpha(i), alpha1(i))
2619      END IF
2620    END DO
2621  END DO
2622
2623  RETURN
2624END SUBROUTINE wake_vec_modulation
2625
2626END MODULE lmdz_wake
Note: See TracBrowser for help on using the repository browser.