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

Last change on this file since 3709 was 3648, checked in by jghattas, 5 years ago

Correction d'erreur introduit en commit [3641] pour faire tourner en mode debug.

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