source: LMDZ6/trunk/libf/phylmd/lmdz_wake.f90 @ 5821

Last change on this file since 5821 was 5804, checked in by fhourdin, 3 months ago

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