source: LMDZ6/branches/Test_modipsl/libf/phylmdiso/wake.F90 @ 5454

Last change on this file since 5454 was 4491, checked in by crisi, 21 months ago

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