source: LMDZ6/trunk/libf/phylmdiso/lmdz_wake.F90 @ 5276

Last change on this file since 5276 was 5276, checked in by abarral, 21 hours ago

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