source: LMDZ6/branches/Amaury_dev/libf/phylmdiso/lmdz_wake.F90 @ 5419

Last change on this file since 5419 was 5160, checked in by abarral, 5 months ago

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