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

Last change on this file since 4283 was 4231, checked in by fhourdin, 21 months ago

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