source: LMDZ6/trunk/libf/phylmdiso/wake.F90 @ 4437

Last change on this file since 4437 was 4374, checked in by lguez, 19 months 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.