source: LMDZ6/branches/Ocean_skin/libf/phylmdiso/wake.F90 @ 4074

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

Initial import of the physics wih isotopes from Camille Risi
CR
File size: 81.7 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
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#ifdef ISOVERIF
500        write(*,*) 'wake 500: debut inits'
501#endif
502!jyg<
503!!  DO k = 1, klev
504!!    DO i = 1, klon
505!!      ppi(i, k) = pi(i, k)
506!!      deltatw0(i, k) = deltatw(i, k)
507!!      deltaqw0(i, k) = deltaqw(i, k)
508!!      te(i, k) = te0(i, k)
509!!      qe(i, k) = qe0(i, k)
510!!      dtls(i, k) = 0.
511!!      dqls(i, k) = 0.
512!!      d_deltat_gw(i, k) = 0.
513!!      d_te(i, k) = 0.
514!!      d_qe(i, k) = 0.
515!!      d_deltatw(i, k) = 0.
516!!      d_deltaqw(i, k) = 0.
517!!      ! IM 060508 beg
518!!      d_deltatw2(i, k) = 0.
519!!      d_deltaqw2(i, k) = 0.
520!!      ! IM 060508 end
521!!    END DO
522!!  END DO
523      ppi(:,:) = pi(:,:)
524      deltatw0(:,:) = deltatw(:,:)
525      deltaqw0(:,:) = deltaqw(:,:)
526      te(:,:) = te0(:,:)
527      qe(:,:) = qe0(:,:)
528      dtls(:,:) = 0.
529      dqls(:,:) = 0.
530      d_deltat_gw(:,:) = 0.
531      d_te(:,:) = 0.
532      d_qe(:,:) = 0.
533      d_deltatw(:,:) = 0.
534      d_deltaqw(:,:) = 0.
535      d_deltatw2(:,:) = 0.
536      d_deltaqw2(:,:) = 0.
537#ifdef ISO
538      deltaxtw0(:,:,:)= deltaxtw(:,:,:)
539      xte(:,:,:) = xte0(:,:,:)
540      dxtls(:,:,:) = 0.
541      d_xte(:,:,:) = 0.
542      d_deltaxtw(:,:,:) = 0.
543      d_deltaxtw2(:,:,:)=0.
544      xt1(:,:,:) = 0.
545      xt2(:,:,:)=0.   
546      ! init non indispensable mais facilite verif iso
547      q1(:,:)=0.0
548      q2(:,:)=0.0
549#endif
550
551#ifdef ISOVERIF
552        write(*,*) 'wake 552: milieu inits'
553#endif
554      IF (iflag_wk_act == 0) THEN
555        act(:) = 0.
556      ELSEIF (iflag_wk_act == 1) THEN
557        act(:) = 1.
558      ENDIF
559
560!!  DO i = 1, klon
561!!   sigmaw_in(i) = sigmaw(i)
562!!  END DO
563   sigmaw_in(:) = sigmaw(:)
564!>jyg
565
566  ! sigmaw1=sigmaw
567  ! IF (sigd_con.GT.sigmaw1) THEN
568  ! print*, 'sigmaw,sigd_con', sigmaw, sigd_con
569  ! ENDIF
570  IF (iflag_wk_pop_dyn >=1) THEN
571    DO i = 1, klon
572      wdens_targ = max(wdens(i),wdensmin)
573      d_wdens2(i) = wdens_targ - wdens(i)
574      wdens(i) = wdens_targ
575    END DO
576  ELSE
577    DO i = 1, klon
578      d_awdens2(i) = 0.
579      d_wdens2(i) = 0.
580    END DO
581  ENDIF  ! (iflag_wk_pop_dyn >=1)
582!
583  DO i = 1, klon
584    ! c      sigmaw(i) = amax1(sigmaw(i),sigd_con(i))
585!jyg<
586!!    sigmaw(i) = amax1(sigmaw(i), sigmad)
587!!    sigmaw(i) = amin1(sigmaw(i), 0.99)
588    sigmaw_targ = min(max(sigmaw(i), sigmad),0.99)
589    d_sigmaw2(i) = sigmaw_targ - sigmaw(i)
590    sigmaw(i) = sigmaw_targ
591!>jyg
592  END DO
593
594!
595  IF (iflag_wk_pop_dyn >= 1) THEN
596    awdens_in(:) = awdens(:)
597    wdens_in(:) = wdens(:)
598!!    wdens(:) = wdens(:) + wgen(:)*dtime
599!!    d_wdens2(:) = wgen(:)*dtime
600!!  ELSE
601  ENDIF  ! (iflag_wk_pop_dyn >= 1)
602
603  wape(:) = 0.
604  wape2(:) = 0.
605  d_sigmaw(:) = 0.
606  ktopw(:) = 0
607!
608!<jyg
609dth(:,:) = 0.
610tu(:,:) = 0.
611qu(:,:) = 0.
612dtke(:,:) = 0.
613dqke(:,:) = 0.
614spread(:,:) = 0.
615omgbdth(:,:) = 0.
616omg(:,:) = 0.
617dp_omgb(:,:) = 0.
618dp_deltomg(:,:) = 0.
619hw(:) = 0.
620wape(:) = 0.
621fip(:) = 0.
622gfl(:) = 0.
623cstar(:) = 0.
624ktopw(:) = 0
625!
626!  Vertical advection local variables
627omgbw(:,:) = 0.
628omgtop(:) = 0
629dp_omgbw(:,:) = 0.
630omgbdq(:,:) = 0.
631#ifdef ISO
632xtu(:,:,:) = 0.
633dxtke(:,:,:) = 0.
634omgbdxt(:,:,:) = 0.
635#endif
636!>jyg
637!
638  IF (prt_level>=10) THEN
639    PRINT *, 'wake-1, sigmaw(igout) ', sigmaw(igout)
640    PRINT *, 'wake-1, deltatw(igout,k) ', (k,deltatw(igout,k), k=1,klev)
641    PRINT *, 'wake-1, deltaqw(igout,k) ', (k,deltaqw(igout,k), k=1,klev)
642    PRINT *, 'wake-1, dowwdraughts, amdwn(igout,k) ', (k,amdwn(igout,k), k=1,klev)
643    PRINT *, 'wake-1, dowwdraughts, dtdwn(igout,k) ', (k,dtdwn(igout,k), k=1,klev)
644    PRINT *, 'wake-1, dowwdraughts, dqdwn(igout,k) ', (k,dqdwn(igout,k), k=1,klev)
645    PRINT *, 'wake-1, updraughts, amup(igout,k) ', (k,amup(igout,k), k=1,klev)
646    PRINT *, 'wake-1, updraughts, dta(igout,k) ', (k,dta(igout,k), k=1,klev)
647    PRINT *, 'wake-1, updraughts, dqa(igout,k) ', (k,dqa(igout,k), k=1,klev)
648  ENDIF
649
650  ! 2. - Prognostic part
651  ! --------------------
652
653
654  ! 2.1 - Undisturbed area and Wake integrals
655  ! ---------------------------------------------------------
656
657  DO i = 1, klon
658    z(i) = 0.
659    ktop(i) = 0
660    kupper(i) = 0
661    sum_thu(i) = 0.
662    sum_tu(i) = 0.
663    sum_qu(i) = 0.
664    sum_thvu(i) = 0.
665    sum_dth(i) = 0.
666    sum_dq(i) = 0.
667    sum_rho(i) = 0.
668    sum_dtdwn(i) = 0.
669    sum_dqdwn(i) = 0.
670
671    av_thu(i) = 0.
672    av_tu(i) = 0.
673    av_qu(i) = 0.
674    av_thvu(i) = 0.
675    av_dth(i) = 0.
676    av_dq(i) = 0.
677    av_rho(i) = 0.
678    av_dtdwn(i) = 0.
679    av_dqdwn(i) = 0.
680    ! pas besoin d'isos ici
681  END DO
682
683
684#ifdef ISOVERIF
685      ! TODO     
686#endif
687
688  ! Distance between wakes
689  DO i = 1, klon
690    ll(i) = (1-sqrt(sigmaw(i)))/sqrt(wdens(i))
691  END DO
692  ! Potential temperatures and humidity
693  ! ----------------------------------------------------------
694  DO k = 1, klev
695    DO i = 1, klon
696      ! write(*,*)'wake 1',i,k,rd,te(i,k)
697      rho(i, k) = p(i, k)/(rd*te(i,k))
698      ! write(*,*)'wake 2',rho(i,k)
699      IF (k==1) THEN
700        ! write(*,*)'wake 3',i,k,rd,te(i,k)
701        rhoh(i, k) = ph(i, k)/(rd*te(i,k))
702        ! write(*,*)'wake 4',i,k,rd,te(i,k)
703        zhh(i, k) = 0
704      ELSE
705        ! write(*,*)'wake 5',rd,(te(i,k)+te(i,k-1))
706        rhoh(i, k) = ph(i, k)*2./(rd*(te(i,k)+te(i,k-1)))
707        ! write(*,*)'wake 6',(-rhoh(i,k)*RG)+zhh(i,k-1)
708        zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*rg) + zhh(i, k-1)
709      END IF
710      ! write(*,*)'wake 7',ppi(i,k)
711      the(i, k) = te(i, k)/ppi(i, k)
712      thu(i, k) = (te(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k)
713      tu(i, k) = te(i, k) - deltatw(i, k)*sigmaw(i)
714      qu(i, k) = qe(i, k) - deltaqw(i, k)*sigmaw(i)
715      ! write(*,*)'wake 8',(rd*(te(i,k)+deltatw(i,k)))
716      rhow(i, k) = p(i, k)/(rd*(te(i,k)+deltatw(i,k)))
717      dth(i, k) = deltatw(i, k)/ppi(i, k)
718#ifdef ISO
719        do ixt=1,ntraciso
720          xtu(ixt,i,k)  =  xte(ixt,i,k) - deltaxtw(ixt,i,k)*sigmaw(i)
721        enddo !do ixt=1,ntraciso
722#endif
723    END DO
724  END DO
725
726
727#ifdef ISO
728#ifdef ISOVERIF
729        ! TODO   
730#endif
731#endif
732
733  DO k = 1, klev - 1
734    DO i = 1, klon
735      IF (k==1) THEN
736        n2(i, k) = 0
737      ELSE
738        n2(i, k) = amax1(0., -rg**2/the(i,k)*rho(i,k)*(the(i,k+1)-the(i,k-1))/ &
739                             (p(i,k+1)-p(i,k-1)))
740      END IF
741      zh(i, k) = (zhh(i,k)+zhh(i,k+1))/2
742
743      cgw(i, k) = sqrt(n2(i,k))*zh(i, k)
744      tgw(i, k) = coefgw*cgw(i, k)/ll(i)
745    END DO
746  END DO
747
748  DO i = 1, klon
749    n2(i, klev) = 0
750    zh(i, klev) = 0
751    cgw(i, klev) = 0
752    tgw(i, klev) = 0
753  END DO
754
755  ! Calcul de la masse volumique moyenne de la colonne   (bdlmd)
756  ! -----------------------------------------------------------------
757
758  DO k = 1, klev
759    DO i = 1, klon
760      epaisseur1(i, k) = 0.
761      epaisseur2(i, k) = 0.
762    END DO
763  END DO
764
765  DO i = 1, klon
766    epaisseur1(i, 1) = -(ph(i,2)-ph(i,1))/(rho(i,1)*rg) + 1.
767    epaisseur2(i, 1) = -(ph(i,2)-ph(i,1))/(rho(i,1)*rg) + 1.
768    rhow_moyen(i, 1) = rhow(i, 1)
769  END DO
770
771  DO k = 2, klev
772    DO i = 1, klon
773      epaisseur1(i, k) = -(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg) + 1.
774      epaisseur2(i, k) = epaisseur2(i, k-1) + epaisseur1(i, k)
775      rhow_moyen(i, k) = (rhow_moyen(i,k-1)*epaisseur2(i,k-1)+rhow(i,k)* &
776        epaisseur1(i,k))/epaisseur2(i, k)
777    END DO
778  END DO
779
780
781  ! Choose an integration bound well above wake top
782  ! -----------------------------------------------------------------
783
784  ! Pupper = 50000.  ! melting level
785  ! Pupper = 60000.
786  ! Pupper = 80000.  ! essais pour case_e
787  DO i = 1, klon
788    pupper(i) = 0.6*ph(i, 1)
789    pupper(i) = max(pupper(i), 45000.)
790    ! cc        Pupper(i) = 60000.
791  END DO
792
793
794  ! Determine Wake top pressure (Ptop) from buoyancy integral
795  ! --------------------------------------------------------
796
797  ! -1/ Pressure of the level where dth becomes less than delta_t_min.
798
799  DO i = 1, klon
800    ptop_provis(i) = ph(i, 1)
801  END DO
802  DO k = 2, klev
803    DO i = 1, klon
804
805      ! IM v3JYG; ptop_provis(i).LT. ph(i,1)
806
807      IF (dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min .AND. &
808          ptop_provis(i)==ph(i,1)) THEN
809        ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)- &
810                          (dth(i,k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
811      END IF
812    END DO
813  END DO
814
815  ! -2/ dth integral
816
817  DO i = 1, klon
818    sum_dth(i) = 0.
819    dthmin(i) = -delta_t_min
820    z(i) = 0.
821  END DO
822
823  DO k = 1, klev
824    DO i = 1, klon
825      dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*rg)
826      IF (dz(i)>0) THEN
827        z(i) = z(i) + dz(i)
828        sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
829        dthmin(i) = amin1(dthmin(i), dth(i,k))
830      END IF
831    END DO
832  END DO
833
834  ! -3/ height of triangle with area= sum_dth and base = dthmin
835
836  DO i = 1, klon
837    hw0(i) = 2.*sum_dth(i)/amin1(dthmin(i), -0.5)
838    hw0(i) = amax1(hwmin, hw0(i))
839  END DO
840
841  ! -4/ now, get Ptop
842
843  DO i = 1, klon
844    z(i) = 0.
845    ptop(i) = ph(i, 1)
846  END DO
847
848  DO k = 1, klev
849    DO i = 1, klon
850      dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg), hw0(i)-z(i))
851      IF (dz(i)>0) THEN
852        z(i) = z(i) + dz(i)
853        ptop(i) = ph(i, k) - rho(i, k)*rg*dz(i)
854      END IF
855    END DO
856  END DO
857
858  IF (prt_level>=10) THEN
859    PRINT *, 'wake-2, ptop_provis(igout), ptop(igout) ', ptop_provis(igout), ptop(igout)
860  ENDIF
861
862
863  ! -5/ Determination de ktop et kupper
864
865  DO k = klev, 1, -1
866    DO i = 1, klon
867      IF (ph(i,k+1)<ptop(i)) ktop(i) = k
868      IF (ph(i,k+1)<pupper(i)) kupper(i) = k
869    END DO
870  END DO
871
872  ! On evite kupper = 1 et kupper = klev
873  DO i = 1, klon
874    kupper(i) = max(kupper(i), 2)
875    kupper(i) = min(kupper(i), klev-1)
876  END DO
877
878
879  ! -6/ Correct ktop and ptop
880
881  DO i = 1, klon
882    ptop_new(i) = ptop(i)
883  END DO
884  DO k = klev, 2, -1
885    DO i = 1, klon
886      IF (k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. &
887          dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
888        ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, &
889          k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
890      END IF
891    END DO
892  END DO
893
894  DO i = 1, klon
895    ptop(i) = ptop_new(i)
896  END DO
897
898  DO k = klev, 1, -1
899    DO i = 1, klon
900      IF (ph(i,k+1)<ptop(i)) ktop(i) = k
901    END DO
902  END DO
903
904  IF (prt_level>=10) THEN
905    PRINT *, 'wake-3, ktop(igout), kupper(igout) ', ktop(igout), kupper(igout)
906  ENDIF
907
908  ! -5/ Set deltatw & deltaqw to 0 above kupper
909
910  DO k = 1, klev
911    DO i = 1, klon
912      IF (k>=kupper(i)) THEN
913        deltatw(i, k) = 0.
914        deltaqw(i, k) = 0.
915        d_deltatw2(i,k) = -deltatw0(i,k)
916#ifdef ISO
917        do ixt=1,ntraciso
918        deltaxtw(ixt,i, k) = 0.
919        d_deltaxtw2(ixt,i,k) = -deltaxtw0(ixt,i,k)
920        enddo
921#endif
922      END IF
923    END DO
924  END DO
925
926
927  ! Vertical gradient of LS omega
928
929  DO k = 1, klev
930    DO i = 1, klon
931      IF (k<=kupper(i)) THEN
932        dp_omgb(i, k) = (omgb(i,k+1)-omgb(i,k))/(ph(i,k+1)-ph(i,k))
933      END IF
934    END DO
935  END DO
936
937  ! Integrals (and wake top level number)
938  ! --------------------------------------
939
940  ! Initialize sum_thvu to 1st level virt. pot. temp.
941
942  DO i = 1, klon
943    z(i) = 1.
944    dz(i) = 1.
945    sum_thvu(i) = thu(i, 1)*(1.+epsim1*qu(i,1))*dz(i)
946    sum_dth(i) = 0.
947  END DO
948
949  DO k = 1, klev
950    DO i = 1, klon
951      dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
952      IF (dz(i)>0) THEN
953        z(i) = z(i) + dz(i)
954        sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i)
955        sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i)
956        sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i)
957        sum_thvu(i) = sum_thvu(i) + thu(i, k)*(1.+epsim1*qu(i,k))*dz(i)
958        sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
959        sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
960        sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i)
961        sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
962        sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
963      END IF
964    END DO
965  END DO
966
967  DO i = 1, klon
968    hw0(i) = z(i)
969  END DO
970
971
972  ! 2.1 - WAPE and mean forcing computation
973  ! ---------------------------------------
974
975  ! ---------------------------------------
976
977  ! Means
978
979  DO i = 1, klon
980    av_thu(i) = sum_thu(i)/hw0(i)
981    av_tu(i) = sum_tu(i)/hw0(i)
982    av_qu(i) = sum_qu(i)/hw0(i)
983    av_thvu(i) = sum_thvu(i)/hw0(i)
984    ! av_thve = sum_thve/hw0
985    av_dth(i) = sum_dth(i)/hw0(i)
986    av_dq(i) = sum_dq(i)/hw0(i)
987    av_rho(i) = sum_rho(i)/hw0(i)
988    av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
989    av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
990
991    wape(i) = -rg*hw0(i)*(av_dth(i)+ &
992        epsim1*(av_thu(i)*av_dq(i)+av_dth(i)*av_qu(i)+av_dth(i)*av_dq(i)))/av_thvu(i)
993
994  END DO
995
996  ! 2.2 Prognostic variable update
997  ! ------------------------------
998
999  ! Filter out bad wakes
1000
1001  DO k = 1, klev
1002    DO i = 1, klon
1003      IF (wape(i)<0.) THEN
1004        deltatw(i, k) = 0.
1005        deltaqw(i, k) = 0.
1006        dth(i, k) = 0.
1007        d_deltatw2(i,k) = -deltatw0(i,k)
1008        d_deltaqw2(i,k) = -deltaqw0(i,k)
1009#ifdef ISO
1010        do ixt=1,ntraciso
1011        deltaxtw(ixt,i, k) = 0.
1012        d_deltaxtw2(ixt,i,k) = -deltaxtw0(ixt,i,k)
1013        enddo
1014#endif
1015      END IF
1016    END DO
1017  END DO
1018
1019  DO i = 1, klon
1020    IF (wape(i)<0.) THEN
1021      wape(i) = 0.
1022      cstar(i) = 0.
1023      hw(i) = hwmin
1024!jyg<
1025!!      sigmaw(i) = amax1(sigmad, sigd_con(i))
1026      sigmaw_targ = max(sigmad, sigd_con(i))
1027      d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
1028      sigmaw(i) = sigmaw_targ
1029!>jyg
1030      fip(i) = 0.
1031      gwake(i) = .FALSE.
1032    ELSE
1033      hw(i) = hw0(i)
1034      cstar(i) = stark*sqrt(2.*wape(i))
1035      gwake(i) = .TRUE.
1036    END IF
1037  END DO
1038
1039
1040  ! Check qx and qw positivity
1041  ! --------------------------
1042  DO i = 1, klon
1043    q0_min(i) = min((qe(i,1)-sigmaw(i)*deltaqw(i,1)),  &
1044                    (qe(i,1)+(1.-sigmaw(i))*deltaqw(i,1)))
1045  END DO
1046  DO k = 2, klev
1047    DO i = 1, klon
1048      q1_min(i) = min((qe(i,k)-sigmaw(i)*deltaqw(i,k)), &
1049                      (qe(i,k)+(1.-sigmaw(i))*deltaqw(i,k)))
1050      IF (q1_min(i)<=q0_min(i)) THEN
1051        q0_min(i) = q1_min(i)
1052      END IF
1053    END DO
1054  END DO
1055
1056  DO i = 1, klon
1057    ok_qx_qw(i) = q0_min(i) >= 0.
1058    alpha(i) = 1.
1059  END DO
1060
1061  IF (prt_level>=10) THEN
1062    PRINT *, 'wake-4, sigmaw(igout), cstar(igout), wape(igout), ktop(igout) ', &
1063                      sigmaw(igout), cstar(igout), wape(igout), ktop(igout)
1064  ENDIF
1065
1066
1067  ! C -----------------------------------------------------------------
1068  ! Sub-time-stepping
1069  ! -----------------
1070
1071  nsub = 10
1072  dtimesub = dtime/nsub
1073
1074  ! ------------------------------------------------------------
1075  DO isubstep = 1, nsub
1076    ! ------------------------------------------------------------
1077
1078    ! wk_adv is the logical flag enabling wake evolution in the time advance
1079    ! loop
1080    DO i = 1, klon
1081      wk_adv(i) = ok_qx_qw(i) .AND. alpha(i) >= 1.
1082    END DO
1083    IF (prt_level>=10) THEN
1084      PRINT *, 'wake-4.1, isubstep,wk_adv(igout),cstar(igout),wape(igout), ptop(igout) ', &
1085                          isubstep,wk_adv(igout),cstar(igout),wape(igout), ptop(igout)
1086    ENDIF
1087
1088    ! cc nrlmd   Ajout d'un recalcul de wdens dans le cas d'un entrainement
1089    ! négatif de ktop à kupper --------
1090    ! cc           On calcule pour cela une densité wdens0 pour laquelle on
1091    ! aurait un entrainement nul ---
1092    !jyg<
1093    ! Dans la configuration avec wdens prognostique, il s'agit d'un cas ou
1094    ! les poches sont insuffisantes pour accueillir tout le flux de masse
1095    ! des descentes unsaturees. Nous faisons alors l'hypothese que la
1096    ! convection profonde cree directement de nouvelles poches, sans passer
1097    ! par les thermiques. La nouvelle valeur de wdens est alors imposée.
1098
1099    DO i = 1, klon
1100      ! c       print *,' isubstep,wk_adv(i),cstar(i),wape(i) ',
1101      ! c     $           isubstep,wk_adv(i),cstar(i),wape(i)
1102      IF (wk_adv(i) .AND. cstar(i)>0.01) THEN
1103        omg(i, kupper(i)+1) = -rg*amdwn(i, kupper(i)+1)/sigmaw(i) + &
1104                               rg*amup(i, kupper(i)+1)/(1.-sigmaw(i))
1105        wdens0 = (sigmaw(i)/(4.*3.14))* &
1106          ((1.-sigmaw(i))*omg(i,kupper(i)+1)/((ph(i,1)-pupper(i))*cstar(i)))**(2)
1107        IF (prt_level >= 10) THEN
1108             print*,'omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i), ph(i,1)-pupper(i)', &
1109                     omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i), ph(i,1)-pupper(i)
1110        ENDIF
1111        IF (wdens(i)<=wdens0*1.1) THEN
1112          IF (iflag_wk_pop_dyn >= 1) THEN
1113             d_wdens2(i) = d_wdens2(i) + wdens0 - wdens(i)
1114          ENDIF
1115          wdens(i) = wdens0
1116        END IF
1117      END IF
1118    END DO
1119
1120    DO i = 1, klon
1121      IF (wk_adv(i)) THEN
1122        gfl(i) = 2.*sqrt(3.14*wdens(i)*sigmaw(i))
1123        rad_wk(i) = sqrt(sigmaw(i)/(3.14*wdens(i)))
1124!jyg<
1125!!        sigmaw(i) = amin1(sigmaw(i), sigmaw_max)
1126        sigmaw_targ = min(sigmaw(i), sigmaw_max)
1127        d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
1128        sigmaw(i) = sigmaw_targ
1129!>jyg
1130      END IF
1131    END DO
1132
1133    IF (iflag_wk_pop_dyn >= 1) THEN
1134!    The variable "death_rate" is significant only when iflag_wk_pop_dyn = 0.
1135!    Here, it has to be set to zero.
1136      death_rate(:) = 0.
1137
1138      IF (iflag_wk_act ==2) THEN
1139      DO i = 1, klon
1140        IF (wk_adv(i)) THEN
1141          wape1_act(i) = abs(cin(i))
1142          wape2_act(i) = 2.*wape1_act(i) + 1.
1143          act(i) = min(1., max(0., (wape(i)-wape1_act(i)) / (wape2_act(i)-wape1_act(i)) ))
1144        ENDIF  ! (wk_adv(i))
1145      ENDDO
1146      ENDIF  ! (iflag_wk_act ==2)
1147
1148
1149      DO i = 1, klon
1150        IF (wk_adv(i)) THEN
1151!!          tau_wk(i) = max(rad_wk(i)/(3.*cstar(i))*((cstar(i)/cstart)**1.5 - 1), 100.)
1152          tau_wk_inv(i) = max( (3.*cstar(i))/(rad_wk(i)*((cstar(i)/cstart)**1.5 - 1)), 0.)
1153          tau_wk_inv_min = min(tau_wk_inv(i), 1./dtimesub)
1154          drdt(i) = (cstar(i) - wgen(i)*(sigmaw(i)/wdens(i)-aa0)/gfl(i)) / &
1155                    (1 + 2*f_shear(i)*(2.*sigmaw(i)-aa0*wdens(i)) - 2.*sigmaw(i))
1156!!                    (1 - 2*sigmaw(i)*(1.-f_shear(i)))
1157          drdt_pos=max(drdt(i),0.)
1158
1159!!          d_wdens(i) = ( wgen(i)*(1.+2.*(sigmaw(i)-sigmad)) &
1160!!                     - wdens(i)*tau_wk_inv_min &
1161!!                     - 2.*gfl(i)*wdens(i)*Cstar(i) )*dtimesub
1162          d_awdens(i) = ( wgen(i) - (1./tau_cv)*(awdens(i) - act(i)*wdens(i)) )*dtimesub
1163          d_wdens(i) = ( wgen(i) - (wdens(i)-awdens(i))*tau_wk_inv_min -  &
1164                         2.*wdens(i)*gfl(i)*drdt_pos )*dtimesub
1165          d_wdens(i) = max(d_wdens(i), wdensmin-wdens(i))
1166
1167!!          d_sigmaw(i) = ( (1.-2*f_shear(i)*sigmaw(i))*(gfl(i)*Cstar(i)+wgen(i)*sigmad/wdens(i)) &
1168!!                      + 2.*f_shear(i)*wgen(i)*sigmaw(i)**2/wdens(i) &
1169!!                      - sigmaw(i)*tau_wk_inv_min )*dtimesub
1170          d_sig_gen(i) = wgen(i)*aa0
1171          d_sig_death(i) = - sigmaw(i)*(1.-awdens(i)/wdens(i))*tau_wk_inv_min
1172!!          d_sig_col(i) = - 2*f_shear(i)*sigmaw(i)*gfl(i)*drdt_pos
1173          d_sig_col(i) = - 2*f_shear(i)*(2.*sigmaw(i)-wdens(i)*aa0)*gfl(i)*drdt_pos
1174          d_sigmaw(i) = ( d_sig_gen(i) + d_sig_death(i) + d_sig_col(i) + gfl(i)*cstar(i) )*dtimesub
1175          d_sigmaw(i) = max(d_sigmaw(i), sigmad-sigmaw(i))
1176        ENDIF
1177      ENDDO
1178
1179      IF (prt_level >= 10) THEN
1180        print *,'wake, cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), drdt(1) ', &
1181                       cstar(1), cstar(1)/cstart, rad_wk(1), tau_wk_inv(1), drdt(1)
1182        print *,'wake, wdens(1), awdens(1), act(1), d_awdens(1) ', &
1183                       wdens(1), awdens(1), act(1), d_awdens(1)
1184        print *,'wake, wgen, -(wdens-awdens)*tau_wk_inv, -2.*wdens*gfl*drdt_pos, d_wdens ', &
1185                       wgen(1), -(wdens(1)-awdens(1))*tau_wk_inv(1), -2.*wdens(1)*gfl(1)*drdt_pos, d_wdens(1)
1186        print *,'wake, d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1) ', &
1187                       d_sig_gen(1), d_sig_death(1), d_sig_col(1), d_sigmaw(1)
1188      ENDIF
1189   
1190    ELSE  !  (iflag_wk_pop_dyn >= 1)
1191
1192    ! cc nrlmd
1193
1194      DO i = 1, klon
1195        IF (wk_adv(i)) THEN
1196          ! cc nrlmd          Introduction du taux de mortalité des poches et
1197          ! test sur sigmaw_max=0.4
1198          ! cc         d_sigmaw(i) = gfl(i)*Cstar(i)*dtimesub
1199          IF (sigmaw(i)>=sigmaw_max) THEN
1200            death_rate(i) = gfl(i)*cstar(i)/sigmaw(i)
1201          ELSE
1202            death_rate(i) = 0.
1203          END IF
1204   
1205          d_sigmaw(i) = gfl(i)*cstar(i)*dtimesub - death_rate(i)*sigmaw(i)* &
1206            dtimesub
1207          ! $              - nat_rate(i)*sigmaw(i)*dtimesub
1208          ! c        print*, 'd_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
1209          ! c     $  death_rate(i),ktop(i),kupper(i)',
1210          ! c     $              d_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
1211          ! c     $  death_rate(i),ktop(i),kupper(i)
1212   
1213          ! sigmaw(i) =sigmaw(i) + gfl(i)*Cstar(i)*dtimesub
1214          ! sigmaw(i) =min(sigmaw(i),0.99)     !!!!!!!!
1215          ! wdens = wdens0/(10.*sigmaw)
1216          ! sigmaw =max(sigmaw,sigd_con)
1217          ! sigmaw =max(sigmaw,sigmad)
1218        END IF
1219      END DO
1220
1221    ENDIF   !  (iflag_wk_pop_dyn >= 1)
1222
1223
1224    ! calcul de la difference de vitesse verticale poche - zone non perturbee
1225    ! IM 060208 differences par rapport au code initial; init. a 0 dp_deltomg
1226    ! IM 060208 et omg sur les niveaux de 1 a klev+1, alors que avant l'on definit
1227    ! IM 060208 au niveau k=1..?
1228    !JYG 161013 Correction : maintenant omg est dimensionne a klev.
1229    DO k = 1, klev
1230      DO i = 1, klon
1231        IF (wk_adv(i)) THEN !!! nrlmd
1232          dp_deltomg(i, k) = 0.
1233        END IF
1234      END DO
1235    END DO
1236    DO k = 1, klev
1237      DO i = 1, klon
1238        IF (wk_adv(i)) THEN !!! nrlmd
1239          omg(i, k) = 0.
1240        END IF
1241      END DO
1242    END DO
1243
1244    DO i = 1, klon
1245      IF (wk_adv(i)) THEN
1246        z(i) = 0.
1247        omg(i, 1) = 0.
1248        dp_deltomg(i, 1) = -(gfl(i)*cstar(i))/(sigmaw(i)*(1-sigmaw(i)))
1249      END IF
1250    END DO
1251
1252    DO k = 2, klev
1253      DO i = 1, klon
1254        IF (wk_adv(i) .AND. k<=ktop(i)) THEN
1255          dz(i) = -(ph(i,k)-ph(i,k-1))/(rho(i,k-1)*rg)
1256          z(i) = z(i) + dz(i)
1257          dp_deltomg(i, k) = dp_deltomg(i, 1)
1258          omg(i, k) = dp_deltomg(i, 1)*z(i)
1259        END IF
1260      END DO
1261    END DO
1262
1263    DO i = 1, klon
1264      IF (wk_adv(i)) THEN
1265        dztop(i) = -(ptop(i)-ph(i,ktop(i)))/(rho(i,ktop(i))*rg)
1266        ztop(i) = z(i) + dztop(i)
1267        omgtop(i) = dp_deltomg(i, 1)*ztop(i)
1268      END IF
1269    END DO
1270
1271    IF (prt_level>=10) THEN
1272      PRINT *, 'wake-4.2, omg(igout,k) ', (k,omg(igout,k), k=1,klev)
1273      PRINT *, 'wake-4.2, omgtop(igout), ptop(igout), ktop(igout) ', &
1274                          omgtop(igout), ptop(igout), ktop(igout)
1275    ENDIF
1276
1277    ! -----------------
1278    ! From m/s to Pa/s
1279    ! -----------------
1280
1281    DO i = 1, klon
1282      IF (wk_adv(i)) THEN
1283        omgtop(i) = -rho(i, ktop(i))*rg*omgtop(i)
1284        dp_deltomg(i, 1) = omgtop(i)/(ptop(i)-ph(i,1))
1285      END IF
1286    END DO
1287
1288    DO k = 1, klev
1289      DO i = 1, klon
1290        IF (wk_adv(i) .AND. k<=ktop(i)) THEN
1291          omg(i, k) = -rho(i, k)*rg*omg(i, k)
1292          dp_deltomg(i, k) = dp_deltomg(i, 1)
1293        END IF
1294      END DO
1295    END DO
1296
1297    ! raccordement lineaire de omg de ptop a pupper
1298
1299    DO i = 1, klon
1300      IF (wk_adv(i) .AND. kupper(i)>ktop(i)) THEN
1301        omg(i, kupper(i)+1) = -rg*amdwn(i, kupper(i)+1)/sigmaw(i) + &
1302          rg*amup(i, kupper(i)+1)/(1.-sigmaw(i))
1303        dp_deltomg(i, kupper(i)) = (omgtop(i)-omg(i,kupper(i)+1))/ &
1304          (ptop(i)-pupper(i))
1305      END IF
1306    END DO
1307
1308    ! c      DO i=1,klon
1309    ! c        print*,'Pente entre 0 et kupper (référence)'
1310    ! c     $           ,omg(i,kupper(i)+1)/(pupper(i)-ph(i,1))
1311    ! c        print*,'Pente entre ktop et kupper'
1312    ! c     $   ,(omg(i,kupper(i)+1)-omgtop(i))/(pupper(i)-ptop(i))
1313    ! c      ENDDO
1314    ! c
1315    DO k = 1, klev
1316      DO i = 1, klon
1317        IF (wk_adv(i) .AND. k>ktop(i) .AND. k<=kupper(i)) THEN
1318          dp_deltomg(i, k) = dp_deltomg(i, kupper(i))
1319          omg(i, k) = omgtop(i) + (ph(i,k)-ptop(i))*dp_deltomg(i, kupper(i))
1320        END IF
1321      END DO
1322    END DO
1323!!    print *,'omg(igout,k) ', (k,omg(igout,k),k=1,klev)
1324    ! cc nrlmd
1325    ! c      DO i=1,klon
1326    ! c      print*,'deltaw_ktop,deltaw_conv',omgtop(i),omg(i,kupper(i)+1)
1327    ! c      END DO
1328    ! cc
1329
1330
1331    ! --    Compute wake average vertical velocity omgbw
1332
1333
1334    DO k = 1, klev
1335      DO i = 1, klon
1336        IF (wk_adv(i)) THEN
1337          omgbw(i, k) = omgb(i, k) + (1.-sigmaw(i))*omg(i, k)
1338        END IF
1339      END DO
1340    END DO
1341    ! --    and its vertical gradient dp_omgbw
1342
1343    DO k = 1, klev-1
1344      DO i = 1, klon
1345        IF (wk_adv(i)) THEN
1346          dp_omgbw(i, k) = (omgbw(i,k+1)-omgbw(i,k))/(ph(i,k+1)-ph(i,k))
1347        END IF
1348      END DO
1349    END DO
1350    DO i = 1, klon
1351      IF (wk_adv(i)) THEN
1352          dp_omgbw(i, klev) = 0.
1353      END IF
1354    END DO
1355
1356    ! --    Upstream coefficients for omgb velocity
1357    ! --    (alpha_up(k) is the coefficient of the value at level k)
1358    ! --    (1-alpha_up(k) is the coefficient of the value at level k-1)
1359    DO k = 1, klev
1360      DO i = 1, klon
1361        IF (wk_adv(i)) THEN
1362          alpha_up(i, k) = 0.
1363          IF (omgb(i,k)>0.) alpha_up(i, k) = 1.
1364        END IF
1365      END DO
1366    END DO
1367
1368    ! Matrix expressing [The,deltatw] from  [Th1,Th2]
1369
1370    DO i = 1, klon
1371      IF (wk_adv(i)) THEN
1372        rre1(i) = 1. - sigmaw(i)
1373        rre2(i) = sigmaw(i)
1374      END IF
1375    END DO
1376    rrd1 = -1.
1377    rrd2 = 1.
1378
1379    ! --    Get [Th1,Th2], dth and [q1,q2]
1380
1381    DO k = 1, klev
1382      DO i = 1, klon
1383        IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN
1384          dth(i, k) = deltatw(i, k)/ppi(i, k)
1385          th1(i, k) = the(i, k) - sigmaw(i)*dth(i, k) ! undisturbed area
1386          th2(i, k) = the(i, k) + (1.-sigmaw(i))*dth(i, k) ! wake
1387          q1(i, k) = qe(i, k) - sigmaw(i)*deltaqw(i, k) ! undisturbed area
1388          q2(i, k) = qe(i, k) + (1.-sigmaw(i))*deltaqw(i, k) ! wake
1389#ifdef ISO
1390          do ixt=1,ntraciso
1391            xt1(ixt,i,k) = xte(ixt,i,k) -sigmaw(i)     *deltaxtw(ixt,i,k) ! undisturbed area
1392            xt2(ixt,i,k) = xte(ixt,i,k) +(1.-sigmaw(i))*deltaxtw(ixt,i,k) ! wake
1393          enddo
1394#endif
1395        END IF
1396      END DO
1397    END DO
1398
1399    DO i = 1, klon
1400      IF (wk_adv(i)) THEN !!! nrlmd
1401        d_th1(i, 1) = 0.
1402        d_th2(i, 1) = 0.
1403        d_dth(i, 1) = 0.
1404        d_q1(i, 1) = 0.
1405        d_q2(i, 1) = 0.
1406        d_dq(i, 1) = 0.
1407#ifdef ISO
1408        do ixt=1,ntraciso
1409          d_xt1(ixt,i,1) = 0.
1410          d_xt2(ixt,i,1) = 0.
1411          d_dxt(ixt,i,1) = 0.
1412        enddo !do ixt=1,ntraciso
1413#endif
1414      END IF
1415    END DO
1416
1417    DO k = 2, klev
1418      DO i = 1, klon
1419        IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN
1420          d_th1(i, k) = th1(i, k-1) - th1(i, k)
1421          d_th2(i, k) = th2(i, k-1) - th2(i, k)
1422          d_dth(i, k) = dth(i, k-1) - dth(i, k)
1423          d_q1(i, k) = q1(i, k-1) - q1(i, k)
1424          d_q2(i, k) = q2(i, k-1) - q2(i, k)
1425          d_dq(i, k) = deltaqw(i, k-1) - deltaqw(i, k)
1426#ifdef ISO
1427        do ixt=1,ntraciso
1428          d_xt1(ixt,i,k) = xt1(ixt,i,k-1)-xt1(ixt,i,k)
1429          d_xt2(ixt,i,k) = xt2(ixt,i,k-1)-xt2(ixt,i,k)
1430          d_dxt(ixt,i,k) = deltaxtw(ixt,i,k-1)-deltaxtw(ixt,i,k)
1431        enddo !do ixt=1,ntraciso
1432#endif
1433        END IF
1434      END DO
1435    END DO
1436
1437    DO i = 1, klon
1438      IF (wk_adv(i)) THEN
1439        omgbdth(i, 1) = 0.
1440        omgbdq(i, 1) = 0.
1441#ifdef ISO
1442        do ixt=1,ntraciso
1443         omgbdxt(ixt,i,1) = 0.
1444        enddo !do ixt=1,ntraciso
1445#endif
1446      END IF
1447    END DO
1448
1449    DO k = 2, klev
1450      DO i = 1, klon
1451        IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN !   loop on interfaces
1452          omgbdth(i, k) = omgb(i, k)*(dth(i,k-1)-dth(i,k))
1453          omgbdq(i, k) = omgb(i, k)*(deltaqw(i,k-1)-deltaqw(i,k))
1454#ifdef ISO
1455        do ixt=1,ntraciso
1456          omgbdxt(ixt,i,k)  = omgb(i,k)*(deltaxtw(ixt,i,k-1) - deltaxtw(ixt,i,k))
1457        enddo !do ixt=1,ntraciso
1458#ifdef ISOVERIF
1459      if (iso_eau.gt.0) then
1460        call iso_verif_egalite(deltaqw(i,k-1),deltaxtw(iso_eau,i,k-1),'wake 1460a')
1461        call iso_verif_egalite(deltaqw(i,k),deltaxtw(iso_eau,i,k),'wake 1460b')
1462        call iso_verif_egalite(omgbdq(i,k),omgbdxt(iso_eau,i,k),'wake 1460c')
1463      endif
1464#endif
1465#endif
1466        END IF
1467      END DO
1468    END DO
1469
1470    IF (prt_level>=10) THEN
1471      PRINT *, 'wake-4.3, th1(igout,k) ', (k,th1(igout,k), k=1,klev)
1472      PRINT *, 'wake-4.3, th2(igout,k) ', (k,th2(igout,k), k=1,klev)
1473      PRINT *, 'wake-4.3, dth(igout,k) ', (k,dth(igout,k), k=1,klev)
1474      PRINT *, 'wake-4.3, omgbdth(igout,k) ', (k,omgbdth(igout,k), k=1,klev)
1475    ENDIF
1476
1477    ! -----------------------------------------------------------------
1478    DO k = 1, klev-1
1479      DO i = 1, klon
1480        IF (wk_adv(i) .AND. k<=kupper(i)-1) THEN
1481          ! -----------------------------------------------------------------
1482
1483          ! Compute redistribution (advective) term
1484
1485          d_deltatw(i, k) = dtimesub/(ph(i,k)-ph(i,k+1))* &
1486            (rrd1*omg(i,k)*sigmaw(i)*d_th1(i,k) - &
1487             rrd2*omg(i,k+1)*(1.-sigmaw(i))*d_th2(i,k+1)- &
1488             (1.-alpha_up(i,k))*omgbdth(i,k)- &
1489             alpha_up(i,k+1)*omgbdth(i,k+1))*ppi(i, k)
1490!           print*,'d_deltatw=', k, d_deltatw(i,k)
1491
1492          d_deltaqw(i, k) = dtimesub/(ph(i,k)-ph(i,k+1))* &
1493            (rrd1*omg(i,k)*sigmaw(i)*d_q1(i,k)- &
1494             rrd2*omg(i,k+1)*(1.-sigmaw(i))*d_q2(i,k+1)- &
1495             (1.-alpha_up(i,k))*omgbdq(i,k)- &
1496             alpha_up(i,k+1)*omgbdq(i,k+1))
1497!           print*,'d_deltaqw=', k, d_deltaqw(i,k)
1498#ifdef ISO
1499        do ixt=1,ntraciso
1500          d_deltaxtw(ixt,i,k) = dtimesub/(Ph(i,k)-Ph(i,k+1))* &
1501            (rrd1*omg(i,k  )*sigmaw(i)     *d_xt1(ixt,i,k)- &
1502            rrd2*omg(i,k+1)*(1.-sigmaw(i))*d_xt2(ixt,i,k+1)- &
1503            (1.-alpha_up(i,k))*omgbdxt(ixt,i,k)- &
1504            alpha_up(i,k+1)*omgbdxt(ixt,i,k+1))
1505        enddo !do ixt=1,ntraciso
1506#ifdef ISOVERIF
1507      if (iso_eau.gt.0) then
1508        call iso_verif_egalite(d_q1(i,k),d_xt1(iso_eau,i,k),'wake 1502a')
1509        call iso_verif_egalite(d_q2(i,k),d_xt2(iso_eau,i,k),'wake 1502b')
1510        call iso_verif_egalite(omgbdq(i,k),omgbdxt(iso_eau,i,k),'wake 1502c')
1511        call iso_verif_egalite(omgbdq(i,k+1),omgbdxt(iso_eau,i,k+1),'wake 1502d')
1512        call iso_verif_egalite(d_deltaqw(i,k),d_deltaxtw(iso_eau,i,k),'wake 1502e')
1513      endif
1514#endif
1515#endif
1516
1517          ! and increment large scale tendencies
1518
1519
1520
1521
1522          ! C
1523          ! -----------------------------------------------------------------
1524          d_te(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_th1(i,k)- &
1525                                  rre2(i)*omg(i,k+1)*(1.-sigmaw(i))*d_th2(i,k+1))/ &
1526                                 (ph(i,k)-ph(i,k+1)) &
1527                                 -sigmaw(i)*(1.-sigmaw(i))*dth(i,k)*(omg(i,k)-omg(i,k+1))/ &
1528                                 (ph(i,k)-ph(i,k+1)) )*ppi(i, k)
1529
1530          d_qe(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_q1(i,k)- &
1531                                  rre2(i)*omg(i,k+1)*(1.-sigmaw(i))*d_q2(i,k+1))/ &
1532                                 (ph(i,k)-ph(i,k+1)) &
1533                                 -sigmaw(i)*(1.-sigmaw(i))*deltaqw(i,k)*(omg(i,k)-omg(i,k+1))/ &
1534                                 (ph(i,k)-ph(i,k+1)) )
1535#ifdef ISO
1536        do ixt=1,ntraciso
1537         d_xte(ixt,i,k) =  dtimesub*( &
1538             ( rre1(i)*omg(i,k  )*sigmaw(i)     *d_xt1(ixt,i,k) &
1539              -rre2(i)*omg(i,k+1)*(1.-sigmaw(i))*d_xt2(ixt,i,k+1) ) &
1540                    /(Ph(i,k)-Ph(i,k+1)) &
1541              -sigmaw(i)*(1.-sigmaw(i))*deltaxtw(ixt,i,k) &
1542              *(omg(i,k)-omg(i,k+1))/(Ph(i,k)-Ph(i,k+1)) &
1543              )
1544        enddo !do ixt=1,ntraciso
1545#endif
1546        ELSE IF (wk_adv(i) .AND. k==kupper(i)) THEN
1547          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)
1548
1549          d_qe(i, k) = dtimesub*(rre1(i)*omg(i,k)*sigmaw(i)*d_q1(i,k)/(ph(i,k)-ph(i,k+1)))
1550
1551#ifdef ISO
1552        do ixt=1,ntraciso
1553         d_xte(ixt,i,k) =  dtimesub*( &
1554            ( rre1(i)*omg(i,k  )*sigmaw(i)  *d_xt1(ixt,i,k)     &   
1555             /(Ph(i,k)-Ph(i,k+1))) &
1556                            )
1557        enddo !do ixt=1,ntraciso
1558#endif
1559        END IF
1560        ! cc
1561      END DO
1562    END DO
1563    ! ------------------------------------------------------------------
1564
1565    IF (prt_level>=10) THEN
1566      PRINT *, 'wake-4.3, d_deltatw(igout,k) ', (k,d_deltatw(igout,k), k=1,klev)
1567      PRINT *, 'wake-4.3, d_deltaqw(igout,k) ', (k,d_deltaqw(igout,k), k=1,klev)
1568    ENDIF
1569
1570    ! Increment state variables
1571!jyg<
1572    IF (iflag_wk_pop_dyn >= 1) THEN
1573      DO k = 1, klev
1574        DO i = 1, klon
1575          IF (wk_adv(i) .AND. k<=kupper(i)) THEN
1576            detr(i,k) = - d_sig_death(i) - d_sig_col(i)     
1577            entr(i,k) = d_sig_gen(i)
1578          ENDIF
1579        ENDDO
1580      ENDDO
1581      ELSE  ! (iflag_wk_pop_dyn >= 1)
1582      DO k = 1, klev
1583        DO i = 1, klon
1584          IF (wk_adv(i) .AND. k<=kupper(i)) THEN
1585            detr(i, k) = 0.
1586   
1587            entr(i, k) = 0.
1588          ENDIF
1589        ENDDO
1590      ENDDO
1591    ENDIF  ! (iflag_wk_pop_dyn >= 1)
1592
1593   
1594
1595    DO k = 1, klev
1596      DO i = 1, klon
1597        ! cc nrlmd       IF( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN
1598        IF (wk_adv(i) .AND. k<=kupper(i)) THEN
1599          ! cc
1600
1601
1602
1603          ! Coefficient de répartition
1604
1605          crep(i, k) = crep_sol*(ph(i,kupper(i))-ph(i,k))/ &
1606            (ph(i,kupper(i))-ph(i,1))
1607          crep(i, k) = crep(i, k) + crep_upper*(ph(i,1)-ph(i,k))/ &
1608            (p(i,1)-ph(i,kupper(i)))
1609
1610
1611          ! Reintroduce compensating subsidence term.
1612
1613          ! dtKE(k)=(dtdwn(k)*Crep(k))/sigmaw
1614          ! dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k))
1615          ! .                   /(1-sigmaw)
1616          ! dqKE(k)=(dqdwn(k)*Crep(k))/sigmaw
1617          ! dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k))
1618          ! .                   /(1-sigmaw)
1619
1620          ! dtKE(k)=(dtdwn(k)*Crep(k)+(1-Crep(k))*dta(k))/sigmaw
1621          ! dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k)*Crep(k))
1622          ! .                   /(1-sigmaw)
1623          ! dqKE(k)=(dqdwn(k)*Crep(k)+(1-Crep(k))*dqa(k))/sigmaw
1624          ! dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k)*Crep(k))
1625          ! .                   /(1-sigmaw)
1626
1627          dtke(i, k) = (dtdwn(i,k)/sigmaw(i)-dta(i,k)/(1.-sigmaw(i)))
1628          dqke(i, k) = (dqdwn(i,k)/sigmaw(i)-dqa(i,k)/(1.-sigmaw(i)))
1629#ifdef ISO
1630        do ixt=1,ntraciso
1631          dxtke(ixt,i,k)=(dxtdwn(ixt,i,k)/sigmaw(i) - dxta(ixt,i,k) &
1632               /(1.-sigmaw(i)))
1633        enddo !do ixt=1,ntraciso
1634#ifdef ISOVERIF
1635      if (iso_eau.gt.0) then
1636        call iso_verif_egalite(dqke(i,k),dxtKE(iso_eau,i,k),'wake 1621a')
1637        call iso_verif_egalite(dqdwn(i,k),dxtdwn(iso_eau,i,k),'wake 1621b')
1638        call iso_verif_egalite(dqa(i,k),dxta(iso_eau,i,k),'wake 1621c')
1639        call iso_verif_egalite(d_deltaqw(i,k),d_deltaxtw(iso_eau,i,k),'wake 1621d')
1640      endif
1641#endif
1642#endif
1643          ! print*,'dtKE= ',dtKE(i,k),' dqKE= ',dqKE(i,k)
1644
1645!
1646
1647          ! cc nrlmd          Prise en compte du taux de mortalité
1648          ! cc               Définitions de entr, detr
1649!jyg<
1650!!            detr(i, k) = 0.
1651!!   
1652!!            entr(i, k) = detr(i, k) + gfl(i)*cstar(i) + &
1653!!              sigmaw(i)*(1.-sigmaw(i))*dp_deltomg(i, k)
1654!!
1655            entr(i, k) = entr(i,k) + gfl(i)*cstar(i) + &
1656                         sigmaw(i)*(1.-sigmaw(i))*dp_deltomg(i, k)   
1657!>jyg
1658            spread(i, k) = (entr(i,k)-detr(i,k))/sigmaw(i)
1659
1660          ! cc        spread(i,k) =
1661          ! (1.-sigmaw(i))*dp_deltomg(i,k)+gfl(i)*Cstar(i)/
1662          ! cc     $  sigmaw(i)
1663
1664
1665          ! ajout d'un effet onde de gravité -Tgw(k)*deltatw(k) 03/02/06 YU
1666          ! Jingmei
1667
1668          ! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltat_gw(i,k),
1669          ! &  Tgw(i,k),deltatw(i,k)
1670          d_deltat_gw(i, k) = d_deltat_gw(i, k) - tgw(i, k)*deltatw(i, k)* &
1671            dtimesub
1672          ! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltatw(i,k)
1673          ff(i) = d_deltatw(i, k)/dtimesub
1674
1675          ! Sans GW
1676
1677          ! deltatw(k)=deltatw(k)+dtimesub*(ff+dtKE(k)-spread(k)*deltatw(k))
1678
1679          ! GW formule 1
1680
1681          ! deltatw(k) = deltatw(k)+dtimesub*
1682          ! $         (ff+dtKE(k) - spread(k)*deltatw(k)-Tgw(k)*deltatw(k))
1683
1684          ! GW formule 2
1685
1686          IF (dtimesub*tgw(i,k)<1.E-10) THEN
1687            d_deltatw(i, k) = dtimesub*(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)) - & ! cc
1690               tgw(i,k)*deltatw(i,k) )
1691          ELSE
1692            d_deltatw(i, k) = 1/tgw(i, k)*(1-exp(-dtimesub*tgw(i,k)))* &
1693               (ff(i)+dtke(i,k) - &
1694                entr(i,k)*deltatw(i,k)/sigmaw(i) - &
1695                (death_rate(i)*sigmaw(i)+detr(i,k))*deltatw(i,k)/(1.-sigmaw(i)) - &
1696                tgw(i,k)*deltatw(i,k) )
1697          END IF
1698
1699          dth(i, k) = deltatw(i, k)/ppi(i, k)
1700
1701          gg(i) = d_deltaqw(i, k)/dtimesub
1702
1703          d_deltaqw(i, k) = dtimesub*(gg(i)+dqke(i,k) - &
1704            entr(i,k)*deltaqw(i,k)/sigmaw(i) - &
1705            (death_rate(i)*sigmaw(i)+detr(i,k))*deltaqw(i,k)/(1.-sigmaw(i)))
1706#ifdef ISO
1707        do ixt=1,ntraciso
1708          gg(i)=d_deltaxtw(ixt,i,k)/dtimesub
1709          d_deltaxtw(ixt,i,k) = dtimesub*(gg(i) + dxtKE(ixt,i,k) - &
1710             entr(i,k)*deltaxtw(ixt,i,k)/sigmaw(i) - &
1711             (death_rate(i)*sigmaw(i)+detr(i,k))*deltaxtw(ixt,i,k)/(1.-sigmaw(i)))
1712        enddo !do ixt=1,ntraciso
1713#ifdef ISOVERIF
1714      if (iso_eau.gt.0) then
1715        call iso_verif_egalite(dqke(i,k),dxtKE(iso_eau,i,k),'wake 1692a')
1716        call iso_verif_egalite(deltaqw(i,k),deltaxtw(iso_eau,i,k),'wake 1692b')
1717        call iso_verif_egalite(d_deltaqw(i,k),d_deltaxtw(iso_eau,i,k),'wake 1692c')
1718      endif
1719#endif
1720#endif
1721          ! cc
1722
1723          ! cc nrlmd
1724          ! cc       d_deltatw2(i,k)=d_deltatw2(i,k)+d_deltatw(i,k)
1725          ! cc       d_deltaqw2(i,k)=d_deltaqw2(i,k)+d_deltaqw(i,k)
1726          ! cc
1727        END IF
1728      END DO
1729    END DO
1730
1731#ifdef ISO
1732#ifdef ISOVERIF
1733      if (iso_eau.gt.0) then
1734             call iso_verif_egalite_vect2D(d_deltaxtw,d_deltaqw, &
1735                 'wake 1359',ntraciso,klon,klev)
1736      endif     
1737#endif         
1738#endif 
1739
1740    ! Scale tendencies so that water vapour remains positive in w and x.
1741
1742    CALL wake_vec_modulation(klon, klev, wk_adv, epsilon, qe, d_qe, deltaqw, &
1743      d_deltaqw, sigmaw, d_sigmaw, alpha)
1744
1745    ! cc nrlmd
1746    ! c      print*,'alpha'
1747    ! c      do i=1,klon
1748    ! c         print*,alpha(i)
1749    ! c      end do
1750    ! cc
1751    DO k = 1, klev
1752      DO i = 1, klon
1753        IF (wk_adv(i) .AND. k<=kupper(i)) THEN
1754          d_te(i, k) = alpha(i)*d_te(i, k)
1755          d_qe(i, k) = alpha(i)*d_qe(i, k)
1756          d_deltatw(i, k) = alpha(i)*d_deltatw(i, k)
1757          d_deltaqw(i, k) = alpha(i)*d_deltaqw(i, k)
1758          d_deltat_gw(i, k) = alpha(i)*d_deltat_gw(i, k)
1759#ifdef ISO
1760        do ixt=1,ntraciso
1761          d_xte(ixt,i,k)=alpha(i)*d_xte(ixt,i,k)
1762          d_deltaxtw(ixt,i,k)=alpha(i)*d_deltaxtw(ixt,i,k)
1763        enddo !do ixt=1,ntraciso
1764#endif
1765        END IF
1766      END DO
1767    END DO
1768    DO i = 1, klon
1769      IF (wk_adv(i)) THEN
1770        d_sigmaw(i) = alpha(i)*d_sigmaw(i)
1771      END IF
1772    END DO
1773
1774    ! Update large scale variables and wake variables
1775    ! IM 060208 manque DO i + remplace DO k=1,kupper(i)
1776    ! IM 060208     DO k = 1,kupper(i)
1777    DO k = 1, klev
1778      DO i = 1, klon
1779        IF (wk_adv(i) .AND. k<=kupper(i)) THEN
1780          dtls(i, k) = dtls(i, k) + d_te(i, k)
1781          dqls(i, k) = dqls(i, k) + d_qe(i, k)
1782#ifdef ISO
1783        do ixt=1,ntraciso
1784          dxtls(ixt,i,k)=dxtls(ixt,i,k)+d_xte(ixt,i,k)
1785        enddo !do ixt=1,ntraciso
1786#endif
1787          ! cc nrlmd
1788          d_deltatw2(i, k) = d_deltatw2(i, k) + d_deltatw(i, k)
1789          d_deltaqw2(i, k) = d_deltaqw2(i, k) + d_deltaqw(i, k)
1790#ifdef ISO
1791        do ixt=1,ntraciso
1792          d_deltaxtw2(ixt,i,k)=d_deltaxtw2(ixt,i,k)+d_deltaxtw(ixt,i,k)
1793        enddo !do ixt=1,ntraciso
1794#endif
1795          ! cc
1796        END IF
1797      END DO
1798    END DO
1799
1800
1801#ifdef ISO
1802#ifdef ISOVERIF
1803      if (iso_eau.gt.0) then
1804        DO k= 1,klev
1805          DO i = 1,klon
1806              call iso_verif_egalite_choix(dxtls(iso_eau,i,k), &
1807                dqls(i,k),'wake 1379',errmax,errmaxrel)
1808          enddo ! DO i = 1,klon     
1809        enddo ! DO k= 1,klev
1810      endif !if (iso_eau.gt.0) then
1811#endif
1812#endif
1813
1814
1815    DO k = 1, klev
1816      DO i = 1, klon
1817        IF (wk_adv(i) .AND. k<=kupper(i)) THEN
1818          te(i, k) = te0(i, k) + dtls(i, k)
1819          qe(i, k) = qe0(i, k) + dqls(i, k)
1820          the(i, k) = te(i, k)/ppi(i, k)
1821          deltatw(i, k) = deltatw(i, k) + d_deltatw(i, k)
1822          deltaqw(i, k) = deltaqw(i, k) + d_deltaqw(i, k)
1823          dth(i, k) = deltatw(i, k)/ppi(i, k)
1824          ! c      print*,'k,qx,qw',k,qe(i,k)-sigmaw(i)*deltaqw(i,k)
1825          ! c     $        ,qe(i,k)+(1-sigmaw(i))*deltaqw(i,k)
1826#ifdef ISO
1827        do ixt=1,ntraciso
1828          xte(ixt,i,k) = xte0(ixt,i,k) + dxtls(ixt,i,k)
1829          deltaxtw(ixt,i,k) = deltaxtw(ixt,i,k)+d_deltaxtw(ixt,i,k)
1830        enddo !do ixt=1,ntraciso
1831#endif
1832        END IF
1833      END DO
1834    END DO
1835!
1836    DO i = 1, klon
1837      IF (wk_adv(i)) THEN
1838        sigmaw(i) = sigmaw(i) + d_sigmaw(i)
1839        d_sigmaw2(i) = d_sigmaw2(i) + d_sigmaw(i)
1840      END IF
1841    END DO
1842
1843#ifdef ISO
1844#ifdef ISOVERIF
1845      if (iso_eau.gt.0) then
1846        DO k= 1,klev
1847          DO i = 1,klon
1848              call iso_verif_egalite_choix(xte(iso_eau,i,k), &
1849                qe(i,k),'wake 1379',errmax,errmaxrel)
1850          enddo ! DO i = 1,klon     
1851        enddo ! DO k= 1,klev
1852      endif !if (iso_eau.gt.0) then     
1853      if (iso_hdo.gt.0) then
1854      call iso_verif_aberrant_enc_vect2D( &
1855            xte,qe, &
1856            'wake 1456, xte apres modifs',ntraciso,klon,klev)
1857!      call iso_verif_aberrant_enc_vect2D_ns(
1858!     :       deltaxtw,deltaqw,
1859!     :       'wake 1518, deltaqw apres modifs',ntraciso,klon,klev)
1860      endif
1861#endif
1862#endif
1863
1864!jyg<
1865    IF (iflag_wk_pop_dyn >= 1) THEN
1866      DO i = 1, klon
1867        IF (wk_adv(i)) THEN
1868          awdens(i) = awdens(i) + d_awdens(i)
1869          wdens(i) = wdens(i) + d_wdens(i)
1870          d_awdens2(i) = d_awdens2(i) + d_awdens(i)
1871          d_wdens2(i) = d_wdens2(i) + d_wdens(i)
1872        END IF
1873      END DO
1874      DO i = 1, klon
1875        IF (wk_adv(i)) THEN
1876          wdens_targ = max(wdens(i),wdensmin)
1877          d_wdens2(i) = d_wdens2(i) + wdens_targ - wdens(i)
1878          wdens(i) = wdens_targ
1879!
1880          wdens_targ = min( max(awdens(i),0.), wdens(i) )
1881          d_awdens2(i) = d_awdens2(i) + wdens_targ - awdens(i)
1882          awdens(i) = wdens_targ
1883        END IF
1884      END DO
1885      DO i = 1, klon
1886        IF (wk_adv(i)) THEN
1887          sigmaw_targ = max(sigmaw(i),sigmad)
1888          d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
1889          sigmaw(i) = sigmaw_targ
1890        END IF
1891      END DO
1892    ENDIF  ! (iflag_wk_pop_dyn >= 1)
1893!>jyg
1894
1895
1896    ! Determine Ptop from buoyancy integral
1897    ! ---------------------------------------
1898
1899    ! -     1/ Pressure of the level where dth changes sign.
1900
1901    DO i = 1, klon
1902      IF (wk_adv(i)) THEN
1903        ptop_provis(i) = ph(i, 1)
1904      END IF
1905    END DO
1906
1907    DO k = 2, klev
1908      DO i = 1, klon
1909        IF (wk_adv(i) .AND. ptop_provis(i)==ph(i,1) .AND. &
1910            dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
1911          ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) - &
1912                            (dth(i,k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
1913        END IF
1914      END DO
1915    END DO
1916
1917    ! -     2/ dth integral
1918
1919    DO i = 1, klon
1920      IF (wk_adv(i)) THEN !!! nrlmd
1921        sum_dth(i) = 0.
1922        dthmin(i) = -delta_t_min
1923        z(i) = 0.
1924      END IF
1925    END DO
1926
1927    DO k = 1, klev
1928      DO i = 1, klon
1929        IF (wk_adv(i)) THEN
1930          dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*rg)
1931          IF (dz(i)>0) THEN
1932            z(i) = z(i) + dz(i)
1933            sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
1934            dthmin(i) = amin1(dthmin(i), dth(i,k))
1935          END IF
1936        END IF
1937      END DO
1938    END DO
1939
1940    ! -     3/ height of triangle with area= sum_dth and base = dthmin
1941
1942    DO i = 1, klon
1943      IF (wk_adv(i)) THEN
1944        hw(i) = 2.*sum_dth(i)/amin1(dthmin(i), -0.5)
1945        hw(i) = amax1(hwmin, hw(i))
1946      END IF
1947    END DO
1948
1949    ! -     4/ now, get Ptop
1950
1951    DO i = 1, klon
1952      IF (wk_adv(i)) THEN !!! nrlmd
1953        ktop(i) = 0
1954        z(i) = 0.
1955      END IF
1956    END DO
1957
1958    DO k = 1, klev
1959      DO i = 1, klon
1960        IF (wk_adv(i)) THEN
1961          dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg), hw(i)-z(i))
1962          IF (dz(i)>0) THEN
1963            z(i) = z(i) + dz(i)
1964            ptop(i) = ph(i, k) - rho(i, k)*rg*dz(i)
1965            ktop(i) = k
1966          END IF
1967        END IF
1968      END DO
1969    END DO
1970
1971    ! 4.5/Correct ktop and ptop
1972
1973    DO i = 1, klon
1974      IF (wk_adv(i)) THEN
1975        ptop_new(i) = ptop(i)
1976      END IF
1977    END DO
1978
1979    DO k = klev, 2, -1
1980      DO i = 1, klon
1981        ! IM v3JYG; IF (k .GE. ktop(i)
1982        IF (wk_adv(i) .AND. k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. &
1983            dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
1984          ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) - &
1985                         (dth(i,k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1))
1986        END IF
1987      END DO
1988    END DO
1989
1990
1991    DO i = 1, klon
1992      IF (wk_adv(i)) THEN
1993        ptop(i) = ptop_new(i)
1994      END IF
1995    END DO
1996
1997    DO k = klev, 1, -1
1998      DO i = 1, klon
1999        IF (wk_adv(i)) THEN !!! nrlmd
2000          IF (ph(i,k+1)<ptop(i)) ktop(i) = k
2001        END IF
2002      END DO
2003    END DO
2004
2005    ! 5/ Set deltatw & deltaqw to 0 above kupper
2006
2007    DO k = 1, klev
2008      DO i = 1, klon
2009        IF (wk_adv(i) .AND. k>=kupper(i)) THEN
2010          deltatw(i, k) = 0.
2011          deltaqw(i, k) = 0.
2012          d_deltatw2(i,k) = -deltatw0(i,k)
2013          d_deltaqw2(i,k) = -deltaqw0(i,k)
2014#ifdef ISO
2015        do ixt=1,ntraciso
2016          deltaxtw(ixt,i,k) = 0.
2017          d_deltaxtw2(ixt,i,k) = -deltaxtw0(ixt,i,k)
2018        enddo !do ixt=1,ntraciso
2019#endif
2020        END IF
2021      END DO
2022    END DO
2023
2024
2025    ! -------------Cstar computation---------------------------------
2026    DO i = 1, klon
2027      IF (wk_adv(i)) THEN !!! nrlmd
2028        sum_thu(i) = 0.
2029        sum_tu(i) = 0.
2030        sum_qu(i) = 0.
2031        sum_thvu(i) = 0.
2032        sum_dth(i) = 0.
2033        sum_dq(i) = 0.
2034        sum_rho(i) = 0.
2035        sum_dtdwn(i) = 0.
2036        sum_dqdwn(i) = 0.
2037
2038        av_thu(i) = 0.
2039        av_tu(i) = 0.
2040        av_qu(i) = 0.
2041        av_thvu(i) = 0.
2042        av_dth(i) = 0.
2043        av_dq(i) = 0.
2044        av_rho(i) = 0.
2045        av_dtdwn(i) = 0.
2046        av_dqdwn(i) = 0.
2047      END IF
2048    END DO
2049
2050    ! Integrals (and wake top level number)
2051    ! --------------------------------------
2052
2053    ! Initialize sum_thvu to 1st level virt. pot. temp.
2054
2055    DO i = 1, klon
2056      IF (wk_adv(i)) THEN !!! nrlmd
2057        z(i) = 1.
2058        dz(i) = 1.
2059        sum_thvu(i) = thu(i, 1)*(1.+epsim1*qu(i,1))*dz(i)
2060        sum_dth(i) = 0.
2061      END IF
2062    END DO
2063
2064    DO k = 1, klev
2065      DO i = 1, klon
2066        IF (wk_adv(i)) THEN !!! nrlmd
2067          dz(i) = -(max(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
2068          IF (dz(i)>0) THEN
2069            z(i) = z(i) + dz(i)
2070            sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i)
2071            sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i)
2072            sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i)
2073            sum_thvu(i) = sum_thvu(i) + thu(i, k)*(1.+epsim1*qu(i,k))*dz(i)
2074            sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
2075            sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
2076            sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i)
2077            sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
2078            sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
2079          END IF
2080        END IF
2081      END DO
2082    END DO
2083
2084    DO i = 1, klon
2085      IF (wk_adv(i)) THEN !!! nrlmd
2086        hw0(i) = z(i)
2087      END IF
2088    END DO
2089
2090
2091    ! - WAPE and mean forcing computation
2092    ! ---------------------------------------
2093
2094    ! ---------------------------------------
2095
2096    ! Means
2097
2098    DO i = 1, klon
2099      IF (wk_adv(i)) THEN !!! nrlmd
2100        av_thu(i) = sum_thu(i)/hw0(i)
2101        av_tu(i) = sum_tu(i)/hw0(i)
2102        av_qu(i) = sum_qu(i)/hw0(i)
2103        av_thvu(i) = sum_thvu(i)/hw0(i)
2104        av_dth(i) = sum_dth(i)/hw0(i)
2105        av_dq(i) = sum_dq(i)/hw0(i)
2106        av_rho(i) = sum_rho(i)/hw0(i)
2107        av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
2108        av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
2109
2110        wape(i) = -rg*hw0(i)*(av_dth(i)+epsim1*(av_thu(i)*av_dq(i) + &
2111                              av_dth(i)*av_qu(i)+av_dth(i)*av_dq(i)))/av_thvu(i)
2112      END IF
2113    END DO
2114
2115    ! Filter out bad wakes
2116
2117    DO k = 1, klev
2118      DO i = 1, klon
2119        IF (wk_adv(i)) THEN !!! nrlmd
2120          IF (wape(i)<0.) THEN
2121            deltatw(i, k) = 0.
2122            deltaqw(i, k) = 0.
2123            dth(i, k) = 0.
2124            d_deltatw2(i,k) = -deltatw0(i,k)
2125            d_deltaqw2(i,k) = -deltaqw0(i,k)
2126#ifdef ISO
2127        do ixt=1,ntraciso
2128          deltaxtw(ixt,i,k) = 0.
2129          d_deltaxtw2(ixt,i,k) = -deltaxtw0(ixt,i,k)
2130        enddo !do ixt=1,ntraciso
2131#endif
2132          END IF
2133        END IF
2134      END DO
2135    END DO
2136
2137    DO i = 1, klon
2138      IF (wk_adv(i)) THEN !!! nrlmd
2139        IF (wape(i)<0.) THEN
2140          wape(i) = 0.
2141          cstar(i) = 0.
2142          hw(i) = hwmin
2143!jyg<
2144!!          sigmaw(i) = max(sigmad, sigd_con(i))
2145          sigmaw_targ = max(sigmad, sigd_con(i))
2146          d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
2147          sigmaw(i) = sigmaw_targ
2148!>jyg
2149          fip(i) = 0.
2150          gwake(i) = .FALSE.
2151        ELSE
2152          cstar(i) = stark*sqrt(2.*wape(i))
2153          gwake(i) = .TRUE.
2154        END IF
2155      END IF
2156    END DO
2157
2158  END DO ! end sub-timestep loop
2159
2160  IF (prt_level>=10) THEN
2161    PRINT *, 'wake-5, sigmaw(igout), cstar(igout), wape(igout), ptop(igout) ', &
2162                      sigmaw(igout), cstar(igout), wape(igout), ptop(igout)
2163  ENDIF
2164
2165
2166  ! ----------------------------------------------------------
2167  ! Determine wake final state; recompute wape, cstar, ktop;
2168  ! filter out bad wakes.
2169  ! ----------------------------------------------------------
2170
2171  ! 2.1 - Undisturbed area and Wake integrals
2172  ! ---------------------------------------------------------
2173
2174  DO i = 1, klon
2175    ! cc nrlmd       if (wk_adv(i)) then !!! nrlmd
2176    IF (ok_qx_qw(i)) THEN
2177      ! cc
2178      z(i) = 0.
2179      sum_thu(i) = 0.
2180      sum_tu(i) = 0.
2181      sum_qu(i) = 0.
2182      sum_thvu(i) = 0.
2183      sum_dth(i) = 0.
2184      sum_half_dth(i) = 0.
2185      sum_dq(i) = 0.
2186      sum_rho(i) = 0.
2187      sum_dtdwn(i) = 0.
2188      sum_dqdwn(i) = 0.
2189
2190      av_thu(i) = 0.
2191      av_tu(i) = 0.
2192      av_qu(i) = 0.
2193      av_thvu(i) = 0.
2194      av_dth(i) = 0.
2195      av_dq(i) = 0.
2196      av_rho(i) = 0.
2197      av_dtdwn(i) = 0.
2198      av_dqdwn(i) = 0.
2199
2200      dthmin(i) = -delta_t_min
2201    END IF
2202  END DO
2203  ! Potential temperatures and humidity
2204  ! ----------------------------------------------------------
2205
2206  DO k = 1, klev
2207    DO i = 1, klon
2208      ! cc nrlmd       IF ( wk_adv(i)) THEN
2209      IF (ok_qx_qw(i)) THEN
2210        ! cc
2211        rho(i, k) = p(i, k)/(rd*te(i,k))
2212        IF (k==1) THEN
2213          rhoh(i, k) = ph(i, k)/(rd*te(i,k))
2214          zhh(i, k) = 0
2215        ELSE
2216          rhoh(i, k) = ph(i, k)*2./(rd*(te(i,k)+te(i,k-1)))
2217          zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*rg) + zhh(i, k-1)
2218        END IF
2219        the(i, k) = te(i, k)/ppi(i, k)
2220        thu(i, k) = (te(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k)
2221        tu(i, k) = te(i, k) - deltatw(i, k)*sigmaw(i)
2222        qu(i, k) = qe(i, k) - deltaqw(i, k)*sigmaw(i)
2223        rhow(i, k) = p(i, k)/(rd*(te(i,k)+deltatw(i,k)))
2224        dth(i, k) = deltatw(i, k)/ppi(i, k)
2225#ifdef ISO
2226        do ixt=1,ntraciso
2227          xtu(ixt,i,k)  =  xte(ixt,i,k) - deltaxtw(ixt,i,k)*sigmaw(i)
2228        enddo !do ixt=1,ntraciso
2229#endif
2230      END IF
2231    END DO
2232  END DO
2233
2234#ifdef ISO
2235#ifdef ISOVERIF
2236      if (iso_hdo.gt.0) then
2237      call iso_verif_aberrant_enc_vect2D( &
2238              xtu,qu, &
2239              'wake 1834, apres modifs',ntraciso,klon,klev)
2240      endif   
2241#endif
2242#endif   
2243
2244  ! Integrals (and wake top level number)
2245  ! -----------------------------------------------------------
2246
2247  ! Initialize sum_thvu to 1st level virt. pot. temp.
2248
2249  DO i = 1, klon
2250    ! cc nrlmd       IF ( wk_adv(i)) THEN
2251    IF (ok_qx_qw(i)) THEN
2252      ! cc
2253      z(i) = 1.
2254      dz(i) = 1.
2255      dz_half(i) = 1.
2256      sum_thvu(i) = thu(i, 1)*(1.+epsim1*qu(i,1))*dz(i)
2257      sum_dth(i) = 0.
2258    END IF
2259  END DO
2260
2261  DO k = 1, klev
2262    DO i = 1, klon
2263      ! cc nrlmd       IF ( wk_adv(i)) THEN
2264      IF (ok_qx_qw(i)) THEN
2265        ! cc
2266        dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
2267        dz_half(i) = -(amax1(ph(i,k+1),0.5*(ptop(i)+ph(i,1)))-ph(i,k))/(rho(i,k)*rg)
2268        IF (dz(i)>0) THEN
2269          z(i) = z(i) + dz(i)
2270          sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i)
2271          sum_tu(i) = sum_tu(i) + tu(i, k)*dz(i)
2272          sum_qu(i) = sum_qu(i) + qu(i, k)*dz(i)
2273          sum_thvu(i) = sum_thvu(i) + thu(i, k)*(1.+epsim1*qu(i,k))*dz(i)
2274          sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i)
2275          sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i)
2276          sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i)
2277          sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i)
2278          sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i)
2279!
2280          dthmin(i) = min(dthmin(i), dth(i,k))
2281        END IF
2282        IF (dz_half(i)>0) THEN
2283          sum_half_dth(i) = sum_half_dth(i) + dth(i, k)*dz_half(i)
2284        END IF
2285      END IF
2286    END DO
2287  END DO
2288
2289  DO i = 1, klon
2290    ! cc nrlmd       IF ( wk_adv(i)) THEN
2291    IF (ok_qx_qw(i)) THEN
2292      ! cc
2293      hw0(i) = z(i)
2294    END IF
2295  END DO
2296
2297  ! - WAPE and mean forcing computation
2298  ! -------------------------------------------------------------
2299
2300  ! Means
2301
2302  DO i = 1, klon
2303    ! cc nrlmd       IF ( wk_adv(i)) THEN
2304    IF (ok_qx_qw(i)) THEN
2305      ! cc
2306      av_thu(i) = sum_thu(i)/hw0(i)
2307      av_tu(i) = sum_tu(i)/hw0(i)
2308      av_qu(i) = sum_qu(i)/hw0(i)
2309      av_thvu(i) = sum_thvu(i)/hw0(i)
2310      av_dth(i) = sum_dth(i)/hw0(i)
2311      av_dq(i) = sum_dq(i)/hw0(i)
2312      av_rho(i) = sum_rho(i)/hw0(i)
2313      av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
2314      av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
2315
2316      wape2(i) = -rg*hw0(i)*(av_dth(i)+epsim1*(av_thu(i)*av_dq(i) + &
2317                             av_dth(i)*av_qu(i)+av_dth(i)*av_dq(i)))/av_thvu(i)
2318    END IF
2319  END DO
2320
2321
2322
2323  ! Prognostic variable update
2324  ! ------------------------------------------------------------
2325
2326  ! Filter out bad wakes
2327
2328  IF (iflag_wk_check_trgl>=1) THEN
2329    ! Check triangular shape of dth profile
2330    DO i = 1, klon
2331      IF (ok_qx_qw(i)) THEN
2332        !! print *,'wake, hw0(i), dthmin(i) ', hw0(i), dthmin(i)
2333        !! print *,'wake, 2.*sum_dth(i)/(hw0(i)*dthmin(i)) ', &
2334        !!                2.*sum_dth(i)/(hw0(i)*dthmin(i))
2335        !! print *,'wake, sum_half_dth(i), sum_dth(i) ', &
2336        !!                sum_half_dth(i), sum_dth(i)
2337        IF ((hw0(i) < 1.) .or. (dthmin(i) >= -delta_t_min) ) THEN
2338          wape2(i) = -1.
2339          !! print *,'wake, rej 1'
2340        ELSE IF (iflag_wk_check_trgl==1.AND.abs(2.*sum_dth(i)/(hw0(i)*dthmin(i)) - 1.) > 0.5) THEN
2341          wape2(i) = -1.
2342          !! print *,'wake, rej 2'
2343        ELSE IF (abs(sum_half_dth(i)) < 0.5*abs(sum_dth(i)) ) THEN
2344          wape2(i) = -1.
2345          !! print *,'wake, rej 3'
2346        END IF
2347      END IF
2348    END DO
2349  END IF
2350
2351
2352  DO k = 1, klev
2353    DO i = 1, klon
2354      ! cc nrlmd        IF ( wk_adv(i) .AND. wape2(i) .LT. 0.) THEN
2355      IF (ok_qx_qw(i) .AND. wape2(i)<0.) THEN
2356        ! cc
2357        deltatw(i, k) = 0.
2358        deltaqw(i, k) = 0.
2359        dth(i, k) = 0.
2360        d_deltatw2(i,k) = -deltatw0(i,k)
2361        d_deltaqw2(i,k) = -deltaqw0(i,k)
2362#ifdef ISO
2363        do ixt=1,ntraciso
2364          deltaxtw(ixt,i,k) = 0.
2365          d_deltaxtw2(ixt,i,k) = -deltaxtw0(ixt,i,k)
2366        enddo !do ixt=1,ntraciso
2367#endif
2368      END IF
2369    END DO
2370  END DO
2371
2372
2373  DO i = 1, klon
2374    ! cc nrlmd       IF ( wk_adv(i)) THEN
2375    IF (ok_qx_qw(i)) THEN
2376      ! cc
2377      IF (wape2(i)<0.) THEN
2378        wape2(i) = 0.
2379        cstar2(i) = 0.
2380        hw(i) = hwmin
2381!jyg<
2382!!      sigmaw(i) = amax1(sigmad, sigd_con(i))
2383      sigmaw_targ = max(sigmad, sigd_con(i))
2384      d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
2385      sigmaw(i) = sigmaw_targ
2386!>jyg
2387        fip(i) = 0.
2388        gwake(i) = .FALSE.
2389      ELSE
2390        IF (prt_level>=10) PRINT *, 'wape2>0'
2391        cstar2(i) = stark*sqrt(2.*wape2(i))
2392        gwake(i) = .TRUE.
2393      END IF
2394    END IF
2395  END DO
2396
2397  DO i = 1, klon
2398    ! cc nrlmd       IF ( wk_adv(i)) THEN
2399    IF (ok_qx_qw(i)) THEN
2400      ! cc
2401      ktopw(i) = ktop(i)
2402    END IF
2403  END DO
2404
2405  DO i = 1, klon
2406    ! cc nrlmd       IF ( wk_adv(i)) THEN
2407    IF (ok_qx_qw(i)) THEN
2408      ! cc
2409      IF (ktopw(i)>0 .AND. gwake(i)) THEN
2410
2411        ! jyg1     Utilisation d'un h_efficace constant ( ~ feeding layer)
2412        ! cc       heff = 600.
2413        ! Utilisation de la hauteur hw
2414        ! c       heff = 0.7*hw
2415        heff(i) = hw(i)
2416
2417        fip(i) = 0.5*rho(i, ktopw(i))*cstar2(i)**3*heff(i)*2* &
2418          sqrt(sigmaw(i)*wdens(i)*3.14)
2419        fip(i) = alpk*fip(i)
2420        ! jyg2
2421      ELSE
2422        fip(i) = 0.
2423      END IF
2424    END IF
2425  END DO
2426
2427  ! Limitation de sigmaw
2428
2429  ! cc nrlmd
2430  ! DO i=1,klon
2431  ! IF (OK_qx_qw(i)) THEN
2432  ! IF (sigmaw(i).GE.sigmaw_max) sigmaw(i)=sigmaw_max
2433  ! ENDIF
2434  ! ENDDO
2435  ! cc
2436
2437  !jyg<
2438  IF (iflag_wk_pop_dyn >= 1) THEN
2439    DO i = 1, klon
2440      kill_wake(i) = ((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
2441          .NOT. ok_qx_qw(i) .OR. (wdens(i) < 2.*wdensmin)
2442    ENDDO
2443  ELSE  ! (iflag_wk_pop_dyn >= 1)
2444    DO i = 1, klon
2445      kill_wake(i) = ((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
2446          .NOT. ok_qx_qw(i)
2447    ENDDO
2448  ENDIF  ! (iflag_wk_pop_dyn >= 1)
2449  !>jyg
2450
2451  DO k = 1, klev
2452    DO i = 1, klon
2453!!jyg      IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
2454!!jyg          .NOT. ok_qx_qw(i)) THEN
2455      IF (kill_wake(i)) THEN
2456        ! cc
2457        dtls(i, k) = 0.
2458        dqls(i, k) = 0.
2459        deltatw(i, k) = 0.
2460        deltaqw(i, k) = 0.
2461        d_deltatw2(i,k) = -deltatw0(i,k)
2462        d_deltaqw2(i,k) = -deltaqw0(i,k)
2463#ifdef ISO
2464        do ixt=1,ntraciso
2465          dxtls(ixt,i,k) = 0.
2466          deltaxtw(ixt,i,k) = 0.
2467          d_deltaxtw2(ixt,i,k) = -deltaxtw0(ixt,i,k)
2468        enddo !do ixt=1,ntraciso
2469#endif
2470      END IF  ! (kill_wake(i))
2471    END DO
2472  END DO
2473
2474  DO i = 1, klon
2475!!jyg    IF (((wape(i)>=wape2(i)) .AND. (wape2(i)<=wapecut)) .OR. (ktopw(i)<=2) .OR. &
2476!!jyg        .NOT. ok_qx_qw(i)) THEN
2477      IF (kill_wake(i)) THEN
2478      ktopw(i) = 0
2479      wape(i) = 0.
2480      cstar(i) = 0.
2481!!jyg   Outside subroutine "Wake" hw, wdens and sigmaw are zero when there are no wakes
2482!!      hw(i) = hwmin                       !jyg
2483!!      sigmaw(i) = sigmad                  !jyg
2484      hw(i) = 0.                            !jyg
2485      fip(i) = 0.
2486!!      sigmaw(i) = 0.                        !jyg
2487      sigmaw_targ = 0.
2488      d_sigmaw2(i) = d_sigmaw2(i) + sigmaw_targ - sigmaw(i)
2489      sigmaw(i) = sigmaw_targ
2490      IF (iflag_wk_pop_dyn >= 1) THEN
2491!!        awdens(i) = 0.
2492!!        wdens(i) = 0.
2493        wdens_targ = 0.
2494        d_wdens2(i) = wdens_targ - wdens(i)
2495        wdens(i) = wdens_targ
2496        wdens_targ = 0.
2497        d_awdens2(i) = wdens_targ - awdens(i)
2498        awdens(i) = wdens_targ
2499      ENDIF  ! (iflag_wk_pop_dyn >= 1)
2500    ELSE  ! (kill_wake(i))
2501      wape(i) = wape2(i)
2502      cstar(i) = cstar2(i)
2503    END IF  ! (kill_wake(i))
2504    ! c        print*,'wape wape2 ktopw OK_qx_qw =',
2505    ! c     $          wape(i),wape2(i),ktopw(i),OK_qx_qw(i)
2506  END DO
2507
2508  IF (prt_level>=10) THEN
2509    PRINT *, 'wake-6, wape wape2 ktopw OK_qx_qw =', &
2510                      wape(igout),wape2(igout),ktopw(igout),OK_qx_qw(igout)
2511  ENDIF
2512
2513
2514  ! -----------------------------------------------------------------
2515  ! Get back to tendencies per second
2516
2517  DO k = 1, klev
2518    DO i = 1, klon
2519
2520      ! cc nrlmd        IF ( wk_adv(i) .AND. k .LE. kupper(i)) THEN
2521!jyg<
2522!!      IF (ok_qx_qw(i) .AND. k<=kupper(i)) THEN
2523      IF (ok_qx_qw(i)) THEN
2524!>jyg
2525        ! cc
2526        dtls(i, k) = dtls(i, k)/dtime
2527        dqls(i, k) = dqls(i, k)/dtime
2528        d_deltatw2(i, k) = d_deltatw2(i, k)/dtime
2529        d_deltaqw2(i, k) = d_deltaqw2(i, k)/dtime
2530        d_deltat_gw(i, k) = d_deltat_gw(i, k)/dtime
2531        ! c      print*,'k,dqls,omg,entr,detr',k,dqls(i,k),omg(i,k),entr(i,k)
2532        ! c     $         ,death_rate(i)*sigmaw(i)
2533#ifdef ISO
2534        do ixt=1,ntraciso
2535          dxtls(ixt,i, k) = dxtls(ixt,i, k)/dtime
2536          d_deltaxtw2(ixt,i, k) = d_deltaxtw2(ixt,i, k)/dtime
2537        enddo
2538#endif
2539      END IF
2540    END DO
2541  END DO
2542!jyg<
2543  DO i = 1, klon
2544    d_sigmaw2(i) = d_sigmaw2(i)/dtime
2545    d_awdens2(i) = d_awdens2(i)/dtime
2546    d_wdens2(i) = d_wdens2(i)/dtime
2547  ENDDO
2548!>jyg
2549
2550
2551
2552  RETURN
2553END SUBROUTINE wake
2554
2555SUBROUTINE wake_vec_modulation(nlon, nl, wk_adv, epsilon, qe, d_qe, deltaqw, &
2556    d_deltaqw, sigmaw, d_sigmaw, alpha)
2557  ! ------------------------------------------------------
2558  ! Dtermination du coefficient alpha tel que les tendances
2559  ! corriges alpha*d_G, pour toutes les grandeurs G, correspondent
2560  ! a une humidite positive dans la zone (x) et dans la zone (w).
2561  ! ------------------------------------------------------
2562  IMPLICIT NONE
2563
2564  ! Input
2565  REAL qe(nlon, nl), d_qe(nlon, nl)
2566  REAL deltaqw(nlon, nl), d_deltaqw(nlon, nl)
2567  REAL sigmaw(nlon), d_sigmaw(nlon)
2568  LOGICAL wk_adv(nlon)
2569  INTEGER nl, nlon
2570  ! Output
2571  REAL alpha(nlon)
2572  ! Internal variables
2573  REAL zeta(nlon, nl)
2574  REAL alpha1(nlon)
2575  REAL x, a, b, c, discrim
2576  REAL epsilon
2577  ! DATA epsilon/1.e-15/
2578  INTEGER i,k
2579
2580  DO k = 1, nl
2581    DO i = 1, nlon
2582      IF (wk_adv(i)) THEN
2583        IF ((deltaqw(i,k)+d_deltaqw(i,k))>=0.) THEN
2584          zeta(i, k) = 0.
2585        ELSE
2586          zeta(i, k) = 1.
2587        END IF
2588      END IF
2589    END DO
2590    DO i = 1, nlon
2591      IF (wk_adv(i)) THEN
2592        x = qe(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + d_qe(i, k) + &
2593          (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - d_sigmaw(i) * &
2594          (deltaqw(i,k)+d_deltaqw(i,k))
2595        a = -d_sigmaw(i)*d_deltaqw(i, k)
2596        b = d_qe(i, k) + (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - &
2597          deltaqw(i, k)*d_sigmaw(i)
2598        c = qe(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + epsilon
2599        discrim = b*b - 4.*a*c
2600        ! print*, 'x, a, b, c, discrim', x, a, b, c, discrim
2601        IF (a+b>=0.) THEN !! Condition suffisante pour la positivité de ovap
2602          alpha1(i) = 1.
2603        ELSE
2604          IF (x>=0.) THEN
2605            alpha1(i) = 1.
2606          ELSE
2607            IF (a>0.) THEN
2608              alpha1(i) = 0.9*min( (2.*c)/(-b+sqrt(discrim)),  &
2609                                   (-b+sqrt(discrim))/(2.*a) )
2610            ELSE IF (a==0.) THEN
2611              alpha1(i) = 0.9*(-c/b)
2612            ELSE
2613              ! print*,'a,b,c discrim',a,b,c discrim
2614              alpha1(i) = 0.9*max( (2.*c)/(-b+sqrt(discrim)),  &
2615                                   (-b+sqrt(discrim))/(2.*a))
2616            END IF
2617          END IF
2618        END IF
2619        alpha(i) = min(alpha(i), alpha1(i))
2620      END IF
2621    END DO
2622  END DO
2623
2624  RETURN
2625END SUBROUTINE wake_vec_modulation
2626
2627
2628
2629
Note: See TracBrowser for help on using the repository browser.