source: LMDZ6/branches/blowing_snow/libf/phylmdiso/wake.F90 @ 5327

Last change on this file since 5327 was 4374, checked in by lguez, 2 years ago

Report modifications from phylmd into phylmdiso

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