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

Last change on this file since 5104 was 5103, checked in by abarral, 4 months ago

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