source: LMDZ6/trunk/libf/phylmdiso/lmdz_wake.F90 @ 5262

Last change on this file since 5262 was 5256, checked in by Laurent Fairhead, 6 weeks ago

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