source: LMDZ6/branches/blowing_snow/libf/phylmd/wake.F90 @ 5378

Last change on this file since 5378 was 4453, checked in by fhourdin, 22 months ago

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