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

Last change on this file since 3288 was 3252, checked in by jyg, 7 years ago

Bug correction concerning the growth rate of the
wake radius.

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