source: LMDZ6/trunk/libf/phylmd/wake.F90 @ 3639

Last change on this file since 3639 was 3455, checked in by jyg, 6 years ago

Adding some forgotten SAVE in wake.F90

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