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

Last change on this file since 4133 was 4089, checked in by fhourdin, 2 years ago

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